From 61acdca09445bbb2eddfab02dfa7dfd9d56d8eb1 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 3 Apr 2026 08:36:47 -0600 Subject: [PATCH 01/26] Added `optika.systems.InterpolatedSystem` as an approximation to `SequentialSystem`. --- optika/systems/_interpolated.py | 73 +++++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) create mode 100644 optika/systems/_interpolated.py diff --git a/optika/systems/_interpolated.py b/optika/systems/_interpolated.py new file mode 100644 index 00000000..fe7acdf1 --- /dev/null +++ b/optika/systems/_interpolated.py @@ -0,0 +1,73 @@ +from typing import Callable, Any +import abc +import dataclasses +import named_arrays as na +import optika +from . import AbstractSystem + +__all__ = [ + "AbstractInterpolatedSystem", +] + + +@dataclasses.dataclass(eq=False, repr=False) +class AbstractInterpolatedSystem( + AbstractSystem, +): + """ + Approximate an exact optical system using interpolation. + """ + + @property + @abc.abstractmethod + def distortion(self) -> Callable[ + [na.SpectralPositionalVectorArray], + na.Cartesian2dVectorArray, + ]: + """ + A distortion model which maps positions on the object plane to + positions on the detector. + """ + + 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 detector plane. + + Parameters + ---------- + coordinates + The vertices of each pixel on the object plane. + axis_wavelength + The logical axis + axis_field + The logical axes corresponding to changing field coordinate. + """ + + position = coordinates.position + + position_new = self.distortion(coordinates) + + result = na.regridding.weights( + coordinates_input=position, + coordinates_output=position_new, + axis_input=axis_field, + axis_output=axis_field, + method="conservative", + ) + + return result + + + def image( + self, + scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], + **kwargs: Any, + ) -> na.SpectralPositionalVectorArray: + pass + From 63d784dd510a8f510316fc9b794b3a0901d49f49 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Wed, 8 Jul 2026 20:26:53 -0600 Subject: [PATCH 02/26] lots of improvements --- optika/systems/__init__.py | 3 + optika/systems/_interpolated.py | 73 -------- optika/systems/_linear.py | 294 ++++++++++++++++++++++++++++++++ 3 files changed, 297 insertions(+), 73 deletions(-) delete mode 100644 optika/systems/_interpolated.py create mode 100644 optika/systems/_linear.py diff --git a/optika/systems/__init__.py b/optika/systems/__init__.py index ce20b15b..ef2a6697 100644 --- a/optika/systems/__init__.py +++ b/optika/systems/__init__.py @@ -4,9 +4,12 @@ from ._systems import AbstractSystem from ._sequential import AbstractSequentialSystem, SequentialSystem +from ._linear import AbstractLinearSystem, LinearSystem __all__ = [ "AbstractSystem", "AbstractSequentialSystem", "SequentialSystem", + "AbstractLinearSystem", + "LinearSystem", ] diff --git a/optika/systems/_interpolated.py b/optika/systems/_interpolated.py deleted file mode 100644 index fe7acdf1..00000000 --- a/optika/systems/_interpolated.py +++ /dev/null @@ -1,73 +0,0 @@ -from typing import Callable, Any -import abc -import dataclasses -import named_arrays as na -import optika -from . import AbstractSystem - -__all__ = [ - "AbstractInterpolatedSystem", -] - - -@dataclasses.dataclass(eq=False, repr=False) -class AbstractInterpolatedSystem( - AbstractSystem, -): - """ - Approximate an exact optical system using interpolation. - """ - - @property - @abc.abstractmethod - def distortion(self) -> Callable[ - [na.SpectralPositionalVectorArray], - na.Cartesian2dVectorArray, - ]: - """ - A distortion model which maps positions on the object plane to - positions on the detector. - """ - - 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 detector plane. - - Parameters - ---------- - coordinates - The vertices of each pixel on the object plane. - axis_wavelength - The logical axis - axis_field - The logical axes corresponding to changing field coordinate. - """ - - position = coordinates.position - - position_new = self.distortion(coordinates) - - result = na.regridding.weights( - coordinates_input=position, - coordinates_output=position_new, - axis_input=axis_field, - axis_output=axis_field, - method="conservative", - ) - - return result - - - def image( - self, - scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], - **kwargs: Any, - ) -> na.SpectralPositionalVectorArray: - pass - diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py new file mode 100644 index 00000000..e023d6c5 --- /dev/null +++ b/optika/systems/_linear.py @@ -0,0 +1,294 @@ +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", +] + + +@dataclasses.dataclass(eq=False, repr=False) +class AbstractLinearSystem( + AbstractSystem, +): + """ + Approximate an exact optical system using a linear forward model. + + Subclasses must define a `distortion` method which maps coordinates on the + object plane to positions on the detector. + Subclasses must also define a `vignetting` method which controls how bright + the scene appears. + """ + + @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.""" + + 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 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.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) + + weights_vignetting = self.vignetting(coordinates_cell) + + weights_stop = self.field_stop( + position=na.Cartesian3dVectorArray( + x=coordinates_cell.position.x, + y=coordinates_cell.position.y, + ), + ) + + weights_area = self.area_effective(coordinates_cell.wavelength) + + weights_input = weights_vignetting * weights_stop * weights_area.value + + result = na.regridding.weights( + coordinates_input=position_sensor, + coordinates_output=self.coordinates_sensor, + axis_input=axis_field, + axis_output=axis_field, + weights_input=weights_input, + method="conservative", + ) + + 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]], + values_input: na.AbstractScalar, + ) -> na.AbstractScalar: + """ + Apply a precomputed set of regridding weights to the photon rate of a + scene, mapping it onto the sensor plane. + + This is the cheap linear part of the forward model: :meth:`weights` + builds the (expensive) conservative-regridding operator once, and this + method reuses it for any scene. + + Parameters + ---------- + weights + The conservative-regridding weights computed by :meth:`weights`. + values_input + The photon rate within each pixel of the object plane. + """ + + values_output = na.regridding.regrid_from_weights( + *weights, + values_input=values_input, + ) + + # restore the unit stripped from `weights_input` inside `weights` + # (the effective collecting area), turning the photon rate per unit + # area into a photon rate per sensor pixel. + values_output = values_output * self.weights_unit + + # photon counts cannot be negative + values_output = np.maximum(values_output, 0) + + return values_output + + def image( + self, + scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + noise: bool = True, + **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. + 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. + noise + Whether to include sensor noise in the result. + kwargs + Additional keyword arguments passed to the sensor's + :meth:`~optika.sensors.AbstractImagingSensor.expose` method, such + as `timedelta`. + """ + + 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: + 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. + volume_wavelength = coordinates.wavelength.volume_cell(axis_wavelength) + volume_field = coordinates.position.volume_cell(axis_field) + volume_field = na.as_named_array(volume_field).cell_centers(axis_wavelength) + volume = volume_wavelength * volume_field + + # energy carried by a single photon at each wavelength + wavelength = coordinates.wavelength.cell_centers(axis_wavelength) + energy_photon = astropy.constants.h * astropy.constants.c / wavelength / u.ph + + # integrate the spectral radiance over each voxel and convert it into + # a photon rate per unit collecting area. `weights` folds in the + # effective area, vignetting, and field stop during regridding. + rate = scene.outputs * volume / energy_photon + rate = rate.to(u.photon / u.s / self.weights_unit) + + weights = self.weights( + coordinates=coordinates, + axis_wavelength=axis_wavelength, + axis_field=axis_field, + ) + + rate_sensor = self.image_from_weights(weights, rate) + + image = na.FunctionArray( + inputs=na.SpectralPositionalVectorArray( + wavelength=coordinates.wavelength, + position=self.coordinates_sensor, + ), + outputs=rate_sensor, + ) + + return self.sensor.expose( + image, + axis_wavelength=axis_wavelength, + noise=noise, + **kwargs, + ) + + +@dataclasses.dataclass(eq=False, repr=False) +class LinearSystem( + AbstractLinearSystem, +): + 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. + """ From ba545351b688e38603bf10af174e3d1c4d8cb86a Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 21 Jul 2026 15:19:30 -0600 Subject: [PATCH 03/26] lots of fixes --- optika/distortion/_distortion.py | 18 ++++++ optika/distortion/_distortion_test.py | 1 + optika/radiometry/_effective_area.py | 9 +++ optika/radiometry/_effective_area_test.py | 1 + optika/radiometry/_vignetting.py | 9 +++ optika/radiometry/_vignetting_test.py | 1 + optika/systems/_linear.py | 72 +++++++++++++++++++---- optika/systems/_sequential.py | 2 + optika/systems/_sequential_test.py | 1 + optika/systems/_systems.py | 2 - optika/systems/_systems_test.py | 1 - 11 files changed, 102 insertions(+), 15 deletions(-) diff --git a/optika/distortion/_distortion.py b/optika/distortion/_distortion.py index 49a7fcaf..0fabb003 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.""" diff --git a/optika/distortion/_distortion_test.py b/optika/distortion/_distortion_test.py index f10185ed..fcb58067 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() 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..beccac2d 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.""" diff --git a/optika/radiometry/_vignetting_test.py b/optika/radiometry/_vignetting_test.py index cbef3679..080cf6da 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() diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index e023d6c5..91cb0723 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -68,6 +68,42 @@ def sensor(self) -> optika.sensors.AbstractImagingSensor: 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`. + + This is the same grid of pixel edges that + :meth:`~optika.sensors.AbstractImagingSensor.collect` bins onto. + """ + sensor = self.sensor + return na.Cartesian2dVectorLinearSpace( + start=sensor.aperture.bound_lower.xy, + stop=sensor.aperture.bound_upper.xy, + axis=sensor.axis_pixel, + num=sensor.num_pixel + 1, + ) + + @property + def shape(self) -> dict[str, int]: + 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, @@ -97,24 +133,34 @@ def weights( # factors at the field cell centers (one fewer point along each axis). coordinates_cell = coordinates.cell_centers(axis_field) - weights_vignetting = self.vignetting(coordinates_cell) - - weights_stop = self.field_stop( - position=na.Cartesian3dVectorArray( - x=coordinates_cell.position.x, - y=coordinates_cell.position.y, - ), - ) + 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.value + 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_field, + axis_output=(axis_pixel.x, axis_pixel.y), weights_input=weights_input, method="conservative", ) @@ -124,7 +170,7 @@ def weights( @property def weights_unit(self): """The units associated with :attr:`weights`.""" - return u.cm ** 2 + return u.cm**2 def image_from_weights( self, @@ -199,6 +245,7 @@ def image( as `timedelta`. """ + scene = scene.explicit coordinates = scene.inputs if axis_field is None: @@ -250,6 +297,7 @@ def image( return self.sensor.expose( image, + direction=self.direction, axis_wavelength=axis_wavelength, noise=noise, **kwargs, @@ -268,7 +316,7 @@ class LinearSystem( 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 @@ -286,7 +334,7 @@ class LinearSystem( 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. diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index a81c8f1c..dce65d70 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -23,6 +23,8 @@ @dataclasses.dataclass(eq=False, repr=False) class AbstractSequentialSystem( + optika.mixins.DxfWritable, + optika.mixins.Plottable, AbstractSystem, ): """ diff --git a/optika/systems/_sequential_test.py b/optika/systems/_sequential_test.py index bd4fa4ce..fee9f6f8 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): diff --git a/optika/systems/_systems.py b/optika/systems/_systems.py index 5a779ea2..14eb4312 100644 --- a/optika/systems/_systems.py +++ b/optika/systems/_systems.py @@ -12,8 +12,6 @@ @dataclasses.dataclass(eq=False, repr=False) class AbstractSystem( - optika.mixins.DxfWritable, - optika.mixins.Plottable, optika.mixins.Printable, optika.mixins.Transformable, optika.mixins.Shaped, 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, From 33b7cd07b9081dc04141731071ccc95be74bcbbe Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 21 Jul 2026 15:19:58 -0600 Subject: [PATCH 04/26] add tests --- optika/systems/_linear_test.py | 135 +++++++++++++++++++++++++++++++++ 1 file changed, 135 insertions(+) create mode 100644 optika/systems/_linear_test.py diff --git a/optika/systems/_linear_test.py b/optika/systems/_linear_test.py new file mode 100644 index 00000000..d8170fbe --- /dev/null +++ b/optika/systems/_linear_test.py @@ -0,0 +1,135 @@ +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=50 * u.arcsec / u.mm, + dispersion=250 * u.nm / u.mm, + angle=0 * u.deg, + reference=na.SpectralPositionalVectorArray( + wavelength=550 * u.nm, + position=na.Cartesian2dVectorArray(0, 0) * u.mm, + ), + ) + + +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() -> 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 * u.W / u.cm**2 / u.arcsec**2 / u.nm, + high=1e-18 * u.W / u.cm**2 / u.arcsec**2 / u.nm, + 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.mm) + + @pytest.mark.parametrize("noise", [False, True]) + def test_image(self, a: optika.systems.AbstractLinearSystem, noise: bool): + scene = _scene() + 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 + + +@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 From d3bf48a73521afd9d9b11ccec9a0a73b6b972f94 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 21 Jul 2026 16:36:41 -0600 Subject: [PATCH 05/26] Document `LinearSystem` with a runnable end-to-end example Rewrite the `AbstractLinearSystem` docstring to describe the full linear forward model (distortion, effective area, and optional vignetting/field stop) instead of the stale "define a distortion/vignetting method" text, and add a `LinearSystem` class docstring with a `jupyter-execute` example that assembles the component models and images a USAF-1951 target, mirroring the `SequentialSystem` example. Also document the `shape` property, clarify that `image()` expects energy spectral radiance, and add a Features bullet for the linear forward model. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- docs/index.rst | 2 + optika/systems/_linear.py | 132 +++++++++++++++++++++++++++++++++++--- 2 files changed, 126 insertions(+), 8 deletions(-) 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/systems/_linear.py b/optika/systems/_linear.py index 91cb0723..beb7a471 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -19,12 +19,19 @@ class AbstractLinearSystem( AbstractSystem, ): """ - Approximate an exact optical system using a linear forward model. - - Subclasses must define a `distortion` method which maps coordinates on the - object plane to positions on the detector. - Subclasses must also define a `vignetting` method which controls how bright - the scene appears. + 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 @@ -87,6 +94,10 @@ def coordinates_sensor(self) -> na.AbstractCartesian2dVectorArray: @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), @@ -225,8 +236,9 @@ def image( Parameters ---------- scene - The spectral radiance of the observed scene, sampled on the - vertices of each pixel on the object plane. + The energy spectral radiance of the observed scene (in units such + as :math:`W / cm^2 / arcsec^2 / nm`), sampled on the vertices of + each pixel on the object plane. axis_wavelength The logical axis of `scene` corresponding to changing wavelength. If :obj:`None` (the default), the single axis of @@ -308,6 +320,110 @@ def image( 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=50 * u.arcsec / u.mm, + dispersion=250 * u.nm / u.mm, + angle=0 * u.deg, + reference=na.SpectralPositionalVectorArray( + wavelength=550 * u.nm, + position=na.Cartesian2dVectorArray(0, 0) * u.mm, + ), + ) + 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.sum("wavelength"), + 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.""" From de8f425cb8995a9437430ad5ebba31a7611f747d Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 21 Jul 2026 20:22:40 -0600 Subject: [PATCH 06/26] Add `SequentialSystem.linearize()` to build a `LinearSystem` Fit the system's distortion, vignetting, and effective area models and assemble them (with the sensor) into a `LinearSystem`. The `pupil` argument is interpreted as cell vertices: the effective area fit uses the vertices to compute pupil cell areas, while the distortion and vignetting fits trace at the corresponding cell centers. The field stop is not modeled, so `field_stop` is left `None`. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_sequential.py | 75 +++++++++++++++++++++++++++++- optika/systems/_sequential_test.py | 46 ++++++++++++++++++ 2 files changed, 120 insertions(+), 1 deletion(-) diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index dce65d70..5987f9f6 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", @@ -1513,6 +1513,79 @@ 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, + ) + return LinearSystem( + area_effective=self.area_effective(pupil=pupil, **kwargs), + distortion=self.distortion(pupil=pupil_centers, degree=degree, **kwargs), + sensor=self.sensor, + vignetting=self.vignetting(pupil=pupil_centers, degree=degree, **kwargs), + ) + def _rayfunction_and_axes( self, wavelength: None | u.Quantity | na.AbstractScalar = None, diff --git a/optika/systems/_sequential_test.py b/optika/systems/_sequential_test.py index fee9f6f8..306278df 100644 --- a/optika/systems/_sequential_test.py +++ b/optika/systems/_sequential_test.py @@ -402,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) From edda0fe5639d207a8f2e4d1479a8228fd8e5b90f Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 21 Jul 2026 20:26:26 -0600 Subject: [PATCH 07/26] fixes --- optika/systems/__init__.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/optika/systems/__init__.py b/optika/systems/__init__.py index ef2a6697..08c58064 100644 --- a/optika/systems/__init__.py +++ b/optika/systems/__init__.py @@ -3,13 +3,13 @@ """ from ._systems import AbstractSystem -from ._sequential import AbstractSequentialSystem, SequentialSystem from ._linear import AbstractLinearSystem, LinearSystem +from ._sequential import AbstractSequentialSystem, SequentialSystem __all__ = [ "AbstractSystem", - "AbstractSequentialSystem", - "SequentialSystem", "AbstractLinearSystem", "LinearSystem", + "AbstractSequentialSystem", + "SequentialSystem", ] From 63810764ca41c5ee56b3eafbd9365123b8b6b619 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 21 Jul 2026 20:52:57 -0600 Subject: [PATCH 08/26] coverage --- optika/systems/_linear.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index beb7a471..28db6cb2 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -265,7 +265,7 @@ def image( if axis_wavelength is None: axis = set(na.shape(coordinates.wavelength)) - set(axis_field) - if len(axis) != 1: + if len(axis) != 1: # pragma: nocover raise ValueError( f"unable to infer `axis_wavelength`: expected exactly one " f"axis of `scene.inputs.wavelength` " From 7cc1aed6c2a5e2622bf6d9cfbd2d97c81f182f54 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Wed, 22 Jul 2026 20:08:18 -0600 Subject: [PATCH 09/26] Add `LinearSystem.backproject` and accept photon or energy radiance `LinearSystem.image` no longer hardcodes an energy-to-photon conversion; it passes the flux through in whatever units the scene radiance is given (energy or photon) and lets the sensor convert, matching `SequentialSystem`. Add `backproject` / `backproject_from_weights`, the transpose of the optical forward model, which projects a detector-plane image back onto the object plane and recovers the scene radiance (relying on the corrected `transpose_weights_conservative`). Fix the `direction` computed by `SequentialSystem.linearize`: use the refracted cosine from `rayfunction_default` (as `sensor.collect` does) rather than the uncalled `rayfunction` method and the raw direction vector, averaging only over the grid axes actually present. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_linear.py | 214 ++++++++++++++++++++++++++++++--- optika/systems/_linear_test.py | 22 +++- optika/systems/_sequential.py | 17 +++ 3 files changed, 233 insertions(+), 20 deletions(-) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index 28db6cb2..630c7ead 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -3,7 +3,6 @@ import dataclasses import numpy as np import astropy.units as u -import astropy.constants import named_arrays as na import optika from . import AbstractSystem @@ -123,7 +122,7 @@ def weights( ) -> 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. + plane to each pixel on the sensor plane. Parameters ---------- @@ -178,6 +177,59 @@ def weights( 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.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.value + + 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`.""" @@ -219,6 +271,40 @@ def image_from_weights( return values_output + def backproject_from_weights( + self, + weights: tuple[na.AbstractScalar, dict[str, int], dict[str, int]], + values_input: na.AbstractScalar, + ) -> na.AbstractScalar: + """ + Apply a precomputed set of transposed regridding weights to a + detector-plane image, spreading it back onto the object plane. + + This is the transpose of :meth:`image_from_weights`: + :meth:`weights_transposed` builds the (expensive) transposed operator + once, and this method reuses it for any image. + + Parameters + ---------- + weights + The transposed regridding weights computed by + :meth:`weights_transposed`. + values_input + The rate within each pixel of the detector plane. + """ + + values_output = na.regridding.regrid_from_weights( + *weights, + values_input=values_input, + ) + + # 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. + values_output = values_output / self.weights_unit + + return values_output + def image( self, scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], @@ -236,9 +322,12 @@ def image( Parameters ---------- scene - The energy spectral radiance of the observed scene (in units such - as :math:`W / cm^2 / arcsec^2 / nm`), sampled on the vertices of - each pixel on the object plane. + 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 @@ -281,15 +370,12 @@ def image( volume_field = na.as_named_array(volume_field).cell_centers(axis_wavelength) volume = volume_wavelength * volume_field - # energy carried by a single photon at each wavelength - wavelength = coordinates.wavelength.cell_centers(axis_wavelength) - energy_photon = astropy.constants.h * astropy.constants.c / wavelength / u.ph - - # integrate the spectral radiance over each voxel and convert it into - # a photon rate per unit collecting area. `weights` folds in the - # effective area, vignetting, and field stop during regridding. - rate = scene.outputs * volume / energy_photon - rate = rate.to(u.photon / u.s / self.weights_unit) + # integrate the spectral radiance over each voxel into a flux per unit + # collecting area. The radiance may be given in either energy or photon + # units; the sensor converts energy to photons if necessary. `weights` + # folds in the effective area, vignetting, and field stop during + # regridding. + rate = scene.outputs * volume weights = self.weights( coordinates=coordinates, @@ -308,13 +394,111 @@ def image( ) return self.sensor.expose( - image, + image=image, direction=self.direction, axis_wavelength=axis_wavelength, noise=noise, **kwargs, ) + def backproject( + self, + image: ( + na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar] + | na.AbstractScalar + ), + coordinates: na.SpectralPositionalVectorArray, + weights: None | tuple[na.AbstractScalar, dict[str, int], dict[str, int]] = None, + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: + """ + Transpose of the linear forward model, :meth:`image`. + + Projects a detector-plane image 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 the optical forward model, the regridding + performed by :meth:`image` before + :meth:`~optika.sensors.AbstractImagingSensor.expose`, and not its + inverse. + + Parameters + ---------- + image + The detector-plane image to project onto the object plane, given as + a rate per sensor pixel (the quantity produced by + :meth:`image_from_weights`). + coordinates + The vertices of each pixel on the object plane to project onto. + 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. + 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. + """ + + 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. + volume_wavelength = coordinates.wavelength.volume_cell(axis_wavelength) + volume_field = coordinates.position.volume_cell(axis_field) + volume_field = na.as_named_array(volume_field).cell_centers(axis_wavelength) + volume = volume_wavelength * volume_field + + 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, + ) + + if isinstance(image, na.AbstractFunctionArray): + image = image.outputs + + radiance = self.backproject_from_weights(weights_transposed, image) + + # recover the spectral radiance by undoing the integration over each + # object-plane voxel performed by :meth:`image`. + radiance = radiance / volume + + return na.FunctionArray( + inputs=coordinates, + outputs=radiance, + ) + @dataclasses.dataclass(eq=False, repr=False) class LinearSystem( diff --git a/optika/systems/_linear_test.py b/optika/systems/_linear_test.py index d8170fbe..6431daa3 100644 --- a/optika/systems/_linear_test.py +++ b/optika/systems/_linear_test.py @@ -54,7 +54,7 @@ def _sensor() -> optika.sensors.ImagingSensor: ) -def _scene() -> na.FunctionArray: +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, @@ -66,8 +66,8 @@ def _scene() -> na.FunctionArray: ), ), outputs=na.random.uniform( - low=0 * u.W / u.cm**2 / u.arcsec**2 / u.nm, - high=1e-18 * u.W / u.cm**2 / u.arcsec**2 / u.nm, + low=0 * radiance, + high=radiance, shape_random=dict(field_x=10, field_y=10), ), ) @@ -93,9 +93,21 @@ def test_coordinates_sensor(self, a: optika.systems.AbstractLinearSystem): assert isinstance(result, na.AbstractCartesian2dVectorArray) assert na.unit_normalized(result).is_equivalent(u.mm) + @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): - scene = _scene() + 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) diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index 5987f9f6..798ec88d 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -1579,10 +1579,27 @@ def linearize( 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), ) From 7befb93db037b62d527cdf43fe64a33b8d34031c Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 10:01:27 -0600 Subject: [PATCH 10/26] Invert `expose` in `LinearSystem.backproject` `backproject` now accepts a detector image of electrons and inverts the sensor response with `ImagingSensor.photons_absorbed` before applying the optical transpose, making `image`/`backproject` a matched pair. Add a `test_backproject` roundtrip. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_linear.py | 30 +++++++++++++++++------------- optika/systems/_linear_test.py | 23 +++++++++++++++++++++++ 2 files changed, 40 insertions(+), 13 deletions(-) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index 630c7ead..21e4e618 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -415,20 +415,19 @@ def backproject( """ Transpose of the linear forward model, :meth:`image`. - Projects a detector-plane image 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 the optical forward model, the regridding - performed by :meth:`image` before - :meth:`~optika.sensors.AbstractImagingSensor.expose`, and not its - inverse. + 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 to project onto the object plane, given as - a rate per sensor pixel (the quantity produced by - :meth:`image_from_weights`). + 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. weights @@ -485,10 +484,15 @@ def backproject( axis_field=axis_field, ) - if isinstance(image, na.AbstractFunctionArray): - image = image.outputs + # invert the detector response, mapping the measured electrons back into + # the photon rate per pixel produced by :meth:`image_from_weights`. + image = self.sensor.photons_absorbed( + image, + direction=self.direction, + axis_wavelength=axis_wavelength, + ) - radiance = self.backproject_from_weights(weights_transposed, image) + radiance = self.backproject_from_weights(weights_transposed, image.outputs) # recover the spectral radiance by undoing the integration over each # object-plane voxel performed by :meth:`image`. diff --git a/optika/systems/_linear_test.py b/optika/systems/_linear_test.py index 6431daa3..7a9e4399 100644 --- a/optika/systems/_linear_test.py +++ b/optika/systems/_linear_test.py @@ -123,6 +123,29 @@ def test_image( assert np.all(result.outputs >= 0 * u.electron) 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) + @pytest.mark.parametrize( argnames="a", From 9d6f0b049866abea6e69a6555b2d11aa3ae1eedc Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 11:00:51 -0600 Subject: [PATCH 11/26] Fold volume and sensor response into `image_from_weights`/`backproject_from_weights` `image_from_weights` now takes a scene and returns electrons (integrating the radiance over each voxel, regridding, and exposing the sensor); `backproject_from_weights` takes the detector image plus object-plane coordinates and returns the backprojected radiance (inverting the sensor response, regridding, and dividing out the voxel volume). `image` and `backproject` become thin wrappers that build the weights and delegate, so callers holding precomputed weights (e.g. ctis) can reuse the full forward/transpose model without recomputing the regridding operator. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_linear.py | 245 ++++++++++++++++++++++++++------------ 1 file changed, 169 insertions(+), 76 deletions(-) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index 21e4e618..2bc9fdec 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -238,72 +238,200 @@ def weights_unit(self): def image_from_weights( self, weights: tuple[na.AbstractScalar, dict[str, int], dict[str, int]], - values_input: na.AbstractScalar, - ) -> na.AbstractScalar: + scene: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + noise: bool = True, + **kwargs: Any, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ - Apply a precomputed set of regridding weights to the photon rate of a - scene, mapping it onto the sensor plane. + Apply a precomputed set of regridding weights to a scene, returning the + electrons measured by the sensor. - This is the cheap linear part of the forward model: :meth:`weights` - builds the (expensive) conservative-regridding operator once, and this - method reuses it for any scene. + 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`. - values_input - The photon rate within each pixel of the object plane. + 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. + noise + Whether to include sensor noise in the result. + kwargs + Additional keyword arguments passed to the sensor's + :meth:`~optika.sensors.AbstractImagingSensor.expose` method, such + as `timedelta`. """ - values_output = na.regridding.regrid_from_weights( + 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. + volume_wavelength = coordinates.wavelength.volume_cell(axis_wavelength) + volume_field = coordinates.position.volume_cell(axis_field) + volume_field = na.as_named_array(volume_field).cell_centers(axis_wavelength) + volume = volume_wavelength * volume_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=values_input, + values_input=rate, ) # restore the unit stripped from `weights_input` inside `weights` - # (the effective collecting area), turning the photon rate per unit - # area into a photon rate per sensor pixel. - values_output = values_output * self.weights_unit + # (the effective collecting area), turning the rate per unit area into + # a rate per sensor pixel. + rate_sensor = rate_sensor * self.weights_unit - # photon counts cannot be negative - values_output = np.maximum(values_output, 0) + # the flux incident on the sensor cannot be negative + rate_sensor = np.maximum(rate_sensor, 0) - return values_output + 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, + **kwargs, + ) def backproject_from_weights( self, weights: tuple[na.AbstractScalar, dict[str, int], dict[str, int]], - values_input: na.AbstractScalar, - ) -> na.AbstractScalar: + image: na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar], + coordinates: na.SpectralPositionalVectorArray, + axis_wavelength: None | str = None, + axis_field: None | tuple[str, str] = None, + ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ Apply a precomputed set of transposed regridding weights to a - detector-plane image, spreading it back onto the object plane. + 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 for any image. + 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`. - values_input - The rate within each pixel of the detector plane. + 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. """ - values_output = na.regridding.regrid_from_weights( + 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. + volume_wavelength = coordinates.wavelength.volume_cell(axis_wavelength) + volume_field = coordinates.position.volume_cell(axis_field) + volume_field = na.as_named_array(volume_field).cell_centers(axis_wavelength) + volume = volume_wavelength * volume_field + + # 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, + ) + + radiance = na.regridding.regrid_from_weights( *weights, - values_input=values_input, + 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. - values_output = values_output / self.weights_unit + radiance = radiance / self.weights_unit - return values_output + # recover the spectral radiance by undoing the integration over each + # object-plane voxel performed by `image`. + radiance = radiance / volume + + return na.FunctionArray( + inputs=coordinates, + outputs=radiance, + ) def image( self, @@ -363,40 +491,20 @@ def image( ) (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. - volume_wavelength = coordinates.wavelength.volume_cell(axis_wavelength) - volume_field = coordinates.position.volume_cell(axis_field) - volume_field = na.as_named_array(volume_field).cell_centers(axis_wavelength) - volume = volume_wavelength * volume_field - - # integrate the spectral radiance over each voxel into a flux per unit - # collecting area. The radiance may be given in either energy or photon - # units; the sensor converts energy to photons if necessary. `weights` - # folds in the effective area, vignetting, and field stop during - # regridding. - rate = scene.outputs * volume - + # `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, ) - rate_sensor = self.image_from_weights(weights, rate) - - 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, + return self.image_from_weights( + weights, + scene, axis_wavelength=axis_wavelength, + axis_field=axis_field, noise=noise, **kwargs, ) @@ -464,13 +572,9 @@ def backproject( ) (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. - volume_wavelength = coordinates.wavelength.volume_cell(axis_wavelength) - volume_field = coordinates.position.volume_cell(axis_field) - volume_field = na.as_named_array(volume_field).cell_centers(axis_wavelength) - volume = volume_wavelength * volume_field - + # `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, @@ -484,23 +588,12 @@ def backproject( axis_field=axis_field, ) - # invert the detector response, mapping the measured electrons back into - # the photon rate per pixel produced by :meth:`image_from_weights`. - image = self.sensor.photons_absorbed( + return self.backproject_from_weights( + weights_transposed, image, - direction=self.direction, + coordinates=coordinates, axis_wavelength=axis_wavelength, - ) - - radiance = self.backproject_from_weights(weights_transposed, image.outputs) - - # recover the spectral radiance by undoing the integration over each - # object-plane voxel performed by :meth:`image`. - radiance = radiance / volume - - return na.FunctionArray( - inputs=coordinates, - outputs=radiance, + axis_field=axis_field, ) From c2f956f5b85471b10ef81b3dc76f19796d303759 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 13:19:56 -0600 Subject: [PATCH 12/26] Add an `uncertainty` flag to `LinearSystem.image` When `uncertainty=True`, `image`/`image_from_weights` attach the standard deviation of the measurement noise to the result as a `NormalUncertainScalarArray`, using the sensor's `uncertainty` method. This lets a caller obtain the expected image and its per-pixel noise width in a single pass, which is needed to compute the uncertainty of a wavelength-integrated image before the integration discards the per-wavelength electrons. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_linear.py | 30 +++++++++++++++++++++++++++++- optika/systems/_linear_test.py | 10 ++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index 2bc9fdec..4c2a248f 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -242,6 +242,7 @@ def image_from_weights( axis_wavelength: None | str = None, axis_field: None | tuple[str, str] = None, noise: bool = True, + uncertainty: bool = False, **kwargs: Any, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ @@ -276,6 +277,11 @@ def image_from_weights( are used. 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 @@ -331,7 +337,7 @@ def image_from_weights( outputs=rate_sensor, ) - return self.sensor.expose( + image = self.sensor.expose( image=image, direction=self.direction, axis_wavelength=axis_wavelength, @@ -339,6 +345,21 @@ def image_from_weights( **kwargs, ) + if uncertainty: + width = self.sensor.uncertainty( + image, + direction=self.direction, + axis_wavelength=axis_wavelength, + ) + image = image.replace( + outputs=na.NormalUncertainScalarArray( + nominal=image.outputs, + width=width.outputs, + ), + ) + + return image + def backproject_from_weights( self, weights: tuple[na.AbstractScalar, dict[str, int], dict[str, int]], @@ -439,6 +460,7 @@ def image( axis_wavelength: None | str = None, axis_field: None | tuple[str, str] = None, noise: bool = True, + uncertainty: bool = False, **kwargs: Any, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ @@ -468,6 +490,11 @@ def image( are used. 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 @@ -506,6 +533,7 @@ def image( axis_wavelength=axis_wavelength, axis_field=axis_field, noise=noise, + uncertainty=uncertainty, **kwargs, ) diff --git a/optika/systems/_linear_test.py b/optika/systems/_linear_test.py index 7a9e4399..48a1b9f1 100644 --- a/optika/systems/_linear_test.py +++ b/optika/systems/_linear_test.py @@ -146,6 +146,16 @@ def test_backproject( assert np.all(np.isfinite(result.outputs.value)) assert result.outputs.sum() > 0 * na.unit(radiance) + 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", From a154d498cf89f740b48197372521368832139707 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 14:23:45 -0600 Subject: [PATCH 13/26] revert to old ruff rules --- pyproject.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 610fdd1e..4b06a09b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -45,3 +45,6 @@ Homepage = "https://github.com/sun-data/optika" Documentation = "https://optika.readthedocs.io/en/latest" [tool.setuptools_scm] + +[lint] +select = ["E4", "E7", "E9", "F"] From 098dbfed690b857f8161f9a93d4bbf4b44710eaa Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 20:23:06 -0600 Subject: [PATCH 14/26] Apply the noise model per readout: add `integrate` to the sensor Read noise is a per-readout effect, but the CTIS `LinearSystem` forward model calls `expose` per wavelength, so adding read noise there over-counted it by a factor of sqrt(N_wavelength). Give `expose`, `uncertainty`, `photons_absorbed`, and `measure` an `integrate` keyword (default `True`): the shot/Fano/QY conversion stays per wavelength (exact, and correct for large EUV quantum yields), then the electrons are summed over wavelength into a single readout and the read noise is applied once. `expose` also gains an `uncertainty` flag so it can return a `NormalUncertainScalarArray` in one pass (evaluating `signal` once). `LinearSystem.image`/`backproject` forward `integrate`; `SequentialSystem.image` now delegates its wavelength integration to `expose` instead of collapsing the wavelength grid itself. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/sensors/_sensors.py | 152 +++++++++++++++++++++++++++----- optika/sensors/_sensors_test.py | 31 ++++++- optika/systems/_linear.py | 44 +++++---- optika/systems/_sequential.py | 10 +-- 4 files changed, 185 insertions(+), 52 deletions(-) diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index 15ba6501..f12ff95b 100644 --- a/optika/sensors/_sensors.py +++ b/optika/sensors/_sensors.py @@ -170,6 +170,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[{axis_wavelength: +0}], + wavelength[{axis_wavelength: ~0}], + ], + 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 +217,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 +234,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 +258,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 +287,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 +335,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 +369,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 +388,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 +411,7 @@ def uncertainty( ], direction: float | na.AbstractScalar = 1, axis_wavelength: None | str = None, + integrate: bool = True, ) -> na.FunctionArray[ na.SpectralPositionalVectorArray, na.AbstractScalar, @@ -330,10 +420,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 +441,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 +462,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 +482,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 +528,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..74ab067e 100644 --- a/optika/sensors/_sensors_test.py +++ b/optika/sensors/_sensors_test.py @@ -121,9 +121,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 +133,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 +174,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/_linear.py b/optika/systems/_linear.py index 4c2a248f..db6ffb88 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -243,6 +243,7 @@ def image_from_weights( axis_field: None | tuple[str, str] = None, noise: bool = True, uncertainty: bool = False, + integrate: bool = True, **kwargs: Any, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ @@ -282,6 +283,11 @@ def image_from_weights( to the result, as a :class:`~named_arrays.NormalUncertainScalarArray`, using the sensor's :meth:`~optika.sensors.AbstractImagingSensor.uncertainty`. + 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`. kwargs Additional keyword arguments passed to the sensor's :meth:`~optika.sensors.AbstractImagingSensor.expose` method, such @@ -337,29 +343,16 @@ def image_from_weights( outputs=rate_sensor, ) - image = self.sensor.expose( + return self.sensor.expose( image=image, direction=self.direction, axis_wavelength=axis_wavelength, noise=noise, + integrate=integrate, + uncertainty=uncertainty, **kwargs, ) - if uncertainty: - width = self.sensor.uncertainty( - image, - direction=self.direction, - axis_wavelength=axis_wavelength, - ) - image = image.replace( - outputs=na.NormalUncertainScalarArray( - nominal=image.outputs, - width=width.outputs, - ), - ) - - return image - def backproject_from_weights( self, weights: tuple[na.AbstractScalar, dict[str, int], dict[str, int]], @@ -367,6 +360,7 @@ def backproject_from_weights( coordinates: na.SpectralPositionalVectorArray, axis_wavelength: None | str = None, axis_field: None | tuple[str, str] = None, + integrate: bool = True, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ Apply a precomputed set of transposed regridding weights to a @@ -402,6 +396,10 @@ def backproject_from_weights( 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`. """ coordinates = coordinates.explicit @@ -433,6 +431,7 @@ def backproject_from_weights( image, direction=self.direction, axis_wavelength=axis_wavelength, + integrate=integrate, ) radiance = na.regridding.regrid_from_weights( @@ -461,6 +460,7 @@ def image( axis_field: None | tuple[str, str] = None, noise: bool = True, uncertainty: bool = False, + integrate: bool = True, **kwargs: Any, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ @@ -495,6 +495,11 @@ def image( to the result, as a :class:`~named_arrays.NormalUncertainScalarArray`, using the sensor's :meth:`~optika.sensors.AbstractImagingSensor.uncertainty`. + 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`. kwargs Additional keyword arguments passed to the sensor's :meth:`~optika.sensors.AbstractImagingSensor.expose` method, such @@ -534,6 +539,7 @@ def image( axis_field=axis_field, noise=noise, uncertainty=uncertainty, + integrate=integrate, **kwargs, ) @@ -547,6 +553,7 @@ def backproject( weights: None | tuple[na.AbstractScalar, dict[str, int], dict[str, int]] = None, axis_wavelength: None | str = None, axis_field: None | tuple[str, str] = None, + integrate: bool = True, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ Transpose of the linear forward model, :meth:`image`. @@ -582,6 +589,10 @@ def backproject( 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`. """ coordinates = coordinates.explicit @@ -622,6 +633,7 @@ def backproject( coordinates=coordinates, axis_wavelength=axis_wavelength, axis_field=axis_field, + integrate=integrate, ) diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index 798ec88d..f43726f3 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -1190,21 +1190,13 @@ 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 distortion( From 455438620eb17f0d838829f5bc5e67aecd243ac1 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 20:45:54 -0600 Subject: [PATCH 15/26] Reconcile `image`/`backproject` signatures across systems Declare the shared forward/transpose interface on `AbstractSystem`: both `image` and `backproject` take the common `scene`/`image`, `axis_wavelength`, `axis_field`, `integrate` (and `noise` for `image`) arguments in the same order. `SequentialSystem.image`, `LinearSystem.image`, and `LinearSystem.image_from_weights`/`backproject_from_weights` are all reordered to that common prefix (in both signatures and docstrings), with their extra arguments (`pupil`/`axis_pupil`, `uncertainty`, `weights`) moved to the end. `SequentialSystem` gains a `backproject` that raises `NotImplementedError`, since a ray-traced system is not a linear operator. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_linear.py | 36 ++++++++++++------------ optika/systems/_sequential.py | 44 +++++++++++++++++++++-------- optika/systems/_systems.py | 52 +++++++++++++++++++++++++++++++++-- 3 files changed, 100 insertions(+), 32 deletions(-) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index db6ffb88..105dc594 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -241,9 +241,9 @@ def image_from_weights( 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, - integrate: bool = True, **kwargs: Any, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ @@ -276,6 +276,11 @@ def image_from_weights( 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 @@ -283,11 +288,6 @@ def image_from_weights( to the result, as a :class:`~named_arrays.NormalUncertainScalarArray`, using the sensor's :meth:`~optika.sensors.AbstractImagingSensor.uncertainty`. - 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`. kwargs Additional keyword arguments passed to the sensor's :meth:`~optika.sensors.AbstractImagingSensor.expose` method, such @@ -458,9 +458,9 @@ def image( 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, - integrate: bool = True, **kwargs: Any, ) -> na.FunctionArray[na.SpectralPositionalVectorArray, na.AbstractScalar]: """ @@ -488,6 +488,11 @@ def image( 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 @@ -495,11 +500,6 @@ def image( to the result, as a :class:`~named_arrays.NormalUncertainScalarArray`, using the sensor's :meth:`~optika.sensors.AbstractImagingSensor.uncertainty`. - 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`. kwargs Additional keyword arguments passed to the sensor's :meth:`~optika.sensors.AbstractImagingSensor.expose` method, such @@ -550,10 +550,10 @@ def backproject( | na.AbstractScalar ), coordinates: na.SpectralPositionalVectorArray, - weights: None | tuple[na.AbstractScalar, dict[str, int], dict[str, int]] = None, axis_wavelength: None | str = None, axis_field: None | tuple[str, str] = None, integrate: bool = True, + 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`. @@ -573,11 +573,6 @@ def backproject( plane, as produced by :meth:`image`. coordinates The vertices of each pixel on the object plane to project onto. - 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. axis_wavelength The logical axis of `coordinates` corresponding to changing wavelength. @@ -593,6 +588,11 @@ def backproject( Whether `image` is a single wavelength-integrated readout, to be spread back over the wavelength bins before the transpose. Defaults to :obj:`True`. + 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 diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index f43726f3..5ac76dbe 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -1090,12 +1090,12 @@ 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. @@ -1107,10 +1107,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`, @@ -1121,12 +1117,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 @@ -1134,6 +1124,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 @@ -1199,6 +1199,26 @@ def image( 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( self, wavelength: None | u.Quantity | na.AbstractScalar = None, diff --git a/optika/systems/_systems.py b/optika/systems/_systems.py index 14eb4312..e53756ea 100644 --- a/optika/systems/_systems.py +++ b/optika/systems/_systems.py @@ -26,17 +26,65 @@ 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. + """ + + @abc.abstractmethod + 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]: + """ + Transpose of the forward model, :meth:`image`. + Maps an image of measured electrons back onto the object plane. + + 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. kwargs Additional keyword arguments used by subclass implementations of this method. From becbe63d4651c18a645d08e2ec925ad2711abf14 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 23 Jul 2026 21:23:24 -0600 Subject: [PATCH 16/26] Express the distortion and sensor grid in pixel coordinates Add `AbstractImagingSensor.pixels()`, mapping an in-plane sensor position to fractional pixel coordinates (pixel 0 at the lower edge of the light-sensitive area). `SequentialSystem.distortion` now fits the sensor-side coordinates in pixels, and `LinearSystem.coordinates_sensor` returns the pixel edges in `u.pix`, so the distortion (plate scale in arcsec/pix, dispersion in nm/pix) and the regrid grid share pixel units. Ray tracing and `collect` keep the physical units they naturally produce. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/sensors/_sensors.py | 21 +++++++++++++++++++++ optika/sensors/_sensors_test.py | 9 +++++++++ optika/systems/_linear.py | 9 +++++---- optika/systems/_linear_test.py | 8 ++++---- optika/systems/_sequential.py | 5 +++-- 5 files changed, 42 insertions(+), 10 deletions(-) diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index f12ff95b..b4265e3e 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 + def collect( self, rays: optika.rays.RayVectorArray, diff --git a/optika/sensors/_sensors_test.py b/optika/sensors/_sensors_test.py index 74ab067e..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=[ diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index 105dc594..9d30dfe8 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -80,13 +80,14 @@ def coordinates_sensor(self) -> na.AbstractCartesian2dVectorArray: The vertices of the sensor pixel grid onto which the scene is regridded, derived from :attr:`sensor`. - This is the same grid of pixel edges that - :meth:`~optika.sensors.AbstractImagingSensor.collect` bins onto. + 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.aperture.bound_lower.xy, - stop=sensor.aperture.bound_upper.xy, + 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, ) diff --git a/optika/systems/_linear_test.py b/optika/systems/_linear_test.py index 48a1b9f1..79630a8d 100644 --- a/optika/systems/_linear_test.py +++ b/optika/systems/_linear_test.py @@ -8,12 +8,12 @@ def _distortion() -> optika.distortion.SimpleDistortionModel: return optika.distortion.SimpleDistortionModel( - plate_scale=50 * u.arcsec / u.mm, - dispersion=250 * u.nm / u.mm, + 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(0, 0) * u.mm, + position=na.Cartesian2dVectorArray(16, 16) * u.pix, ), ) @@ -91,7 +91,7 @@ def test_sensor(self, a: optika.systems.AbstractLinearSystem): 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.mm) + assert na.unit_normalized(result).is_equivalent(u.pix) @pytest.mark.parametrize( argnames="radiance", diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index 5ac76dbe..acaa316e 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -1282,9 +1282,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, ) From 257688ff8080fb839a448b410fe6a4a4290c8f16 Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Fri, 24 Jul 2026 07:30:50 -0600 Subject: [PATCH 17/26] Fit the polynomial distortion and vignetting models per channel `PolynomialDistortionModel.fit`/`fit_inverse` and `PolynomialVignettingModel.fit` did not pass `axis_polynomial` to `named_arrays.PolynomialFitFunctionArray.from_degree`, so the least-squares sums ran over every axis of the calibration points, including axes orthogonal to the scene such as the channel axis of a multi-channel instrument. All channels were collapsed into a single averaged polynomial that fit none of them. Passing the scene axes restricts the fit sums to wavelength/field, giving an independent polynomial per orthogonal element. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01HuYXL1BkWsWGAxYnpD8Kqz --- optika/distortion/_distortion.py | 2 ++ optika/distortion/_distortion_test.py | 40 +++++++++++++++++++++++++++ optika/radiometry/_vignetting.py | 1 + optika/radiometry/_vignetting_test.py | 25 +++++++++++++++++ 4 files changed, 68 insertions(+) diff --git a/optika/distortion/_distortion.py b/optika/distortion/_distortion.py index 0fabb003..3116c4e9 100644 --- a/optika/distortion/_distortion.py +++ b/optika/distortion/_distortion.py @@ -376,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, ) @@ -392,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 fcb58067..182656d2 100644 --- a/optika/distortion/_distortion_test.py +++ b/optika/distortion/_distortion_test.py @@ -152,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/_vignetting.py b/optika/radiometry/_vignetting.py index beccac2d..2098ce2a 100644 --- a/optika/radiometry/_vignetting.py +++ b/optika/radiometry/_vignetting.py @@ -184,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 080cf6da..43f7a9db 100644 --- a/optika/radiometry/_vignetting_test.py +++ b/optika/radiometry/_vignetting_test.py @@ -126,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) From c4731c83635c34625c4111e9779691c65d5d1f94 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 24 Jul 2026 19:19:16 -0600 Subject: [PATCH 18/26] Fix docs build and input handling for the pixel-coordinate LinearSystem - Update the LinearSystem and SequentialSystem docstring examples for the pixel-based sensor grid: express the SimpleDistortionModel in pix instead of mm (so it matches the pixel `coordinates_sensor`) and plot the image directly, since the `integrate=True` default already collapses the wavelength axis (dropping the now-invalid `.sum("wavelength")`). - Coerce `coordinates` to a plain SpectralPositionalVectorArray in `weights()`/`weights_transposed()` so richer scene vectors (e.g. Doppler scenes) are normalized to wavelength+position before distortion. - Normalize `AbstractImagingSensor.pixels()` output to pix units. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/sensors/_sensors.py | 2 +- optika/systems/_linear.py | 18 ++++++++++++++---- optika/systems/_sequential.py | 2 +- 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index b4265e3e..aa27105f 100644 --- a/optika/sensors/_sensors.py +++ b/optika/sensors/_sensors.py @@ -108,7 +108,7 @@ def pixels( units. """ lower = self.aperture.bound_lower.xy - return (position - lower) / self.width_pixel * u.pix + return (position - lower) / self.width_pixel * u.pix << u.pix def collect( self, diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index 9d30dfe8..db217a9c 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -135,6 +135,11 @@ def weights( The logical axes corresponding to changing field coordinate. """ + coordinates = na.SpectralPositionalVectorArray( + wavelength=coordinates.wavelength, + position=coordinates.position, + ) + coordinates = coordinates.cell_centers(axis_wavelength) position_sensor = self.distortion.distort(coordinates).position @@ -199,6 +204,11 @@ def weights_transposed( The logical axes corresponding to changing field coordinate. """ + coordinates = na.SpectralPositionalVectorArray( + wavelength=coordinates.wavelength, + position=coordinates.position, + ) + coordinates = coordinates.cell_centers(axis_wavelength) position_sensor = self.distortion.distort(coordinates).position @@ -671,12 +681,12 @@ class LinearSystem( # The distortion, effective area, and sensor models that # define the system. distortion = optika.distortion.SimpleDistortionModel( - plate_scale=50 * u.arcsec / u.mm, - dispersion=250 * u.nm / u.mm, + 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(0, 0) * u.mm, + position=na.Cartesian2dVectorArray(32, 32) * u.pix, ), ) area_effective = optika.radiometry.InterpolatedEffectiveAreaModel( @@ -739,7 +749,7 @@ class LinearSystem( ) na.plt.pcolormesh( image.inputs.position, - C=image.outputs.value.sum("wavelength"), + C=image.outputs.value, ax=ax[1], ) ax[0].set_title("scene") diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index acaa316e..e548231b 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -2117,7 +2117,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( From b3a074f5e63f22e4d8ee7fe80ba0077c08de7a35 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 24 Jul 2026 19:30:53 -0600 Subject: [PATCH 19/26] Cover axis inference in image_from_weights/backproject_from_weights Codecov flagged 6 uncovered lines in optika/systems/_linear.py: the `axis_field`/`axis_wavelength` inference branches of `image_from_weights` and `backproject_from_weights`, which never ran because `image()`/ `backproject()` always pass those axes explicitly. Add `test_from_weights_infer_axes`, which calls the `*_from_weights` methods directly with the axes left to infer from the scene, bringing the module to 100%. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_linear_test.py | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/optika/systems/_linear_test.py b/optika/systems/_linear_test.py index 79630a8d..efccf5a8 100644 --- a/optika/systems/_linear_test.py +++ b/optika/systems/_linear_test.py @@ -146,6 +146,32 @@ def test_backproject( assert np.all(np.isfinite(result.outputs.value)) assert result.outputs.sum() > 0 * na.unit(radiance) + 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_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) From 3b75e6ef830c52a05c263d56f6c069f2ced489ed Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Mon, 27 Jul 2026 14:50:02 -0600 Subject: [PATCH 20/26] Let `backproject` express the radiance in photon or energy units The forward model, `image`, accepts a scene in either photon or energy units (the sensor detects and converts), so make the transpose symmetric: add a `unit` parameter to `AbstractSystem.backproject`, `LinearSystem.backproject`, and `backproject_from_weights` so the caller can request the backprojected spectral radiance in whichever units match the scene. The conversion (`_radiance_to_unit`) scales the natural photon radiance by the energy per photon (hc/lambda) when an energy unit is requested, using the wavelength already available in `backproject_from_weights`. The default (`unit=None`) leaves the radiance in photon units, unchanged from before. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_linear.py | 53 ++++++++++++++++++++++++++++++++++ optika/systems/_linear_test.py | 19 ++++++++++++ optika/systems/_systems.py | 7 +++++ 3 files changed, 79 insertions(+) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index db217a9c..10d8af49 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -3,6 +3,7 @@ import dataclasses import numpy as np import astropy.units as u +import astropy.constants import named_arrays as na import optika from . import AbstractSystem @@ -13,6 +14,31 @@ ] +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, @@ -372,6 +398,7 @@ def backproject_from_weights( 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 @@ -411,6 +438,14 @@ def backproject_from_weights( 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 @@ -459,6 +494,14 @@ def backproject_from_weights( # 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, @@ -564,6 +607,7 @@ def backproject( 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]: """ @@ -599,6 +643,14 @@ def backproject( 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`. @@ -645,6 +697,7 @@ def backproject( axis_wavelength=axis_wavelength, axis_field=axis_field, integrate=integrate, + unit=unit, ) diff --git a/optika/systems/_linear_test.py b/optika/systems/_linear_test.py index efccf5a8..84393371 100644 --- a/optika/systems/_linear_test.py +++ b/optika/systems/_linear_test.py @@ -146,6 +146,25 @@ def test_backproject( assert np.all(np.isfinite(result.outputs.value)) assert result.outputs.sum() > 0 * na.unit(radiance) + 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. diff --git a/optika/systems/_systems.py b/optika/systems/_systems.py index e53756ea..fb98289e 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 @@ -64,6 +65,7 @@ def backproject( 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]: """ @@ -85,6 +87,11 @@ def backproject( 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. From 45ec42d57ce24fbaf2cbe27887838fe1874cb1b0 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Tue, 28 Jul 2026 09:40:09 -0600 Subject: [PATCH 21/26] Use `SpectralPositionalVectorArray.volume_cell` for the voxel volume Replace the hand-written `wavelength.volume_cell(...) * position.volume_cell(...)` (with the manual cell-center alignment) in `image_from_weights` and `backproject_from_weights` with the new `coordinates.volume_cell(...)` from named-arrays 2.3. Bump the named-arrays lower bound to ~=2.3 accordingly. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_linear.py | 10 ++-------- pyproject.toml | 2 +- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index 10d8af49..02f9a6b1 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -350,10 +350,7 @@ def image_from_weights( # volume of each voxel on the object plane: the spectral bin width # times the solid angle (or area) subtended by each field pixel. - volume_wavelength = coordinates.wavelength.volume_cell(axis_wavelength) - volume_field = coordinates.position.volume_cell(axis_field) - volume_field = na.as_named_array(volume_field).cell_centers(axis_wavelength) - volume = volume_wavelength * volume_field + volume = coordinates.volume_cell((axis_wavelength, *axis_field)) # integrate the spectral radiance over each voxel into a flux per unit # collecting area. @@ -466,10 +463,7 @@ def backproject_from_weights( # volume of each voxel on the object plane: the spectral bin width # times the solid angle (or area) subtended by each field pixel. - volume_wavelength = coordinates.wavelength.volume_cell(axis_wavelength) - volume_field = coordinates.position.volume_cell(axis_field) - volume_field = na.as_named_array(volume_field).cell_centers(axis_wavelength) - volume = volume_wavelength * volume_field + volume = coordinates.volume_cell((axis_wavelength, *axis_field)) # invert the detector response, mapping the measured electrons back into # the photon rate per pixel produced by `image_from_weights`. diff --git a/pyproject.toml b/pyproject.toml index 66005d0e..80cff5d2 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,7 +15,7 @@ classifiers = [ ] dependencies = [ "astropy!=6.1.5", - "named-arrays~=2.1", + "named-arrays~=2.3", "pymupdf", "joblib", "ezdxf", From 5d624a56223066041b2fa1e731b00ed0f5285c4e Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 30 Jul 2026 11:23:30 -0600 Subject: [PATCH 22/26] Accept a Doppler object-plane grid in image and backproject `image_from_weights` and `backproject_from_weights` computed the voxel volume with `coordinates.volume_cell(...)`, which raised `NotImplementedError` for a `DopplerPositionalVectorArray` scene (as produced by the ctis `IdealInstrument`) since `volume_cell` is only defined for the spectral-positional vector. Coerce to a plain `SpectralPositionalVectorArray` just for the volume computation, mirroring the ctis `_volume_scene` helper. The coercion is local to the volume, so `image` is unchanged for a normal scene and `backproject` returns the radiance on the original grid (a Doppler grid stays a Doppler grid, preserving `wavelength_rest`). Add `test_image_doppler` and `test_backproject_doppler` (with a `_scene_doppler` helper) covering the Doppler path across both linear system configurations. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_linear.py | 20 ++++++++++++---- optika/systems/_linear_test.py | 44 ++++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+), 4 deletions(-) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index 02f9a6b1..1f0f15b3 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -349,8 +349,14 @@ def image_from_weights( (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. - volume = coordinates.volume_cell((axis_wavelength, *axis_field)) + # times the solid angle (or area) subtended by each field pixel. Coerce + # to a plain spectral-positional vector so a scene defined on a Doppler + # grid (as produced by the ctis `IdealInstrument`) is accepted; only the + # observed wavelength and position of each voxel are needed. + volume = na.SpectralPositionalVectorArray( + wavelength=coordinates.wavelength, + position=coordinates.position, + ).volume_cell((axis_wavelength, *axis_field)) # integrate the spectral radiance over each voxel into a flux per unit # collecting area. @@ -462,8 +468,14 @@ def backproject_from_weights( (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. - volume = coordinates.volume_cell((axis_wavelength, *axis_field)) + # times the solid angle (or area) subtended by each field pixel. Coerce + # to a plain spectral-positional vector so a Doppler grid (as used by + # the ctis `IdealInstrument`) is accepted; the radiance is returned on + # the original grid. + volume = na.SpectralPositionalVectorArray( + wavelength=coordinates.wavelength, + position=coordinates.position, + ).volume_cell((axis_wavelength, *axis_field)) # invert the detector response, mapping the measured electrons back into # the photon rate per pixel produced by `image_from_weights`. diff --git a/optika/systems/_linear_test.py b/optika/systems/_linear_test.py index 84393371..6277f22e 100644 --- a/optika/systems/_linear_test.py +++ b/optika/systems/_linear_test.py @@ -73,6 +73,26 @@ def _scene(radiance: u.Quantity | na.AbstractScalar) -> na.FunctionArray: ) +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, ): @@ -123,6 +143,17 @@ def test_image( 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=[ @@ -146,6 +177,19 @@ def test_backproject( 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) From f7014060628f2df5e0ca41b01f5f3930f39e7b91 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 31 Jul 2026 11:09:15 -0600 Subject: [PATCH 23/26] Use `.spectral_positional` to accept Doppler grids Replace the hand-written `na.SpectralPositionalVectorArray(wavelength=..., position=...)` coercions in `LinearSystem` (`weights`, `weights_transposed`, `image_from_weights`, `backproject_from_weights`) with the `.spectral_positional` property added in named-arrays 2.4.0, and bump the pin accordingly. Behavior is unchanged; the property centralizes the spectral-/Doppler-positional normalization. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/systems/_linear.py | 38 ++++++++++++++------------------------ pyproject.toml | 2 +- 2 files changed, 15 insertions(+), 25 deletions(-) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index 1f0f15b3..c07606e3 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -161,10 +161,7 @@ def weights( The logical axes corresponding to changing field coordinate. """ - coordinates = na.SpectralPositionalVectorArray( - wavelength=coordinates.wavelength, - position=coordinates.position, - ) + coordinates = coordinates.spectral_positional coordinates = coordinates.cell_centers(axis_wavelength) @@ -230,10 +227,7 @@ def weights_transposed( The logical axes corresponding to changing field coordinate. """ - coordinates = na.SpectralPositionalVectorArray( - wavelength=coordinates.wavelength, - position=coordinates.position, - ) + coordinates = coordinates.spectral_positional coordinates = coordinates.cell_centers(axis_wavelength) @@ -349,14 +343,12 @@ def image_from_weights( (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. Coerce - # to a plain spectral-positional vector so a scene defined on a Doppler - # grid (as produced by the ctis `IdealInstrument`) is accepted; only the - # observed wavelength and position of each voxel are needed. - volume = na.SpectralPositionalVectorArray( - wavelength=coordinates.wavelength, - position=coordinates.position, - ).volume_cell((axis_wavelength, *axis_field)) + # 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. @@ -468,14 +460,12 @@ def backproject_from_weights( (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. Coerce - # to a plain spectral-positional vector so a Doppler grid (as used by - # the ctis `IdealInstrument`) is accepted; the radiance is returned on - # the original grid. - volume = na.SpectralPositionalVectorArray( - wavelength=coordinates.wavelength, - position=coordinates.position, - ).volume_cell((axis_wavelength, *axis_field)) + # 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) + ) # invert the detector response, mapping the measured electrons back into # the photon rate per pixel produced by `image_from_weights`. diff --git a/pyproject.toml b/pyproject.toml index 80cff5d2..63cc8109 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -15,7 +15,7 @@ classifiers = [ ] dependencies = [ "astropy!=6.1.5", - "named-arrays~=2.3", + "named-arrays~=2.4", "pymupdf", "joblib", "ezdxf", From 532e288524e329ac261561e04a848f8f80e0570d Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Fri, 31 Jul 2026 12:50:45 -0600 Subject: [PATCH 24/26] Convert the effective area to `weights_unit` before stripping its unit (#194) `weights` and `weights_transposed` stripped the unit of the effective area with `.value`, while `image_from_weights` and `backproject_from_weights` restore it as the hardcoded `weights_unit` (cm^2). If the effective area model evaluates in any other unit (mm^2 for the ESIS-II model), every image is wrong by the unit ratio (x100). Convert with `.to_value(self.weights_unit)` instead, and add a regression test checking that the image is invariant under a unit change of the effective area model. Co-authored-by: Claude Fable 5 --- optika/systems/_linear.py | 6 ++++-- optika/systems/_linear_test.py | 21 +++++++++++++++++++++ 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index c07606e3..6e2dccba 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -191,7 +191,9 @@ def weights( weights_area = self.area_effective(coordinates_cell.wavelength) - weights_input = weights_vignetting * weights_stop * weights_area.value + weights_input = ( + weights_vignetting * weights_stop * weights_area.to_value(self.weights_unit) + ) axis_pixel = self.sensor.axis_pixel @@ -246,7 +248,7 @@ def weights_transposed( weights_area = self.area_effective(coordinates_cell.wavelength) - weights_input = weights_vignetting * weights_area.value + weights_input = weights_vignetting * weights_area.to_value(self.weights_unit) axis_pixel = self.sensor.axis_pixel diff --git a/optika/systems/_linear_test.py b/optika/systems/_linear_test.py index 6277f22e..c35d8b9a 100644 --- a/optika/systems/_linear_test.py +++ b/optika/systems/_linear_test.py @@ -1,3 +1,5 @@ +import dataclasses + import pytest import numpy as np import astropy.units as u @@ -235,6 +237,25 @@ def test_from_weights_infer_axes(self, a: optika.systems.AbstractLinearSystem): ) 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) From 2ace99a5640b5326bbeb62a522398eae5eada45e Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 31 Jul 2026 13:15:53 -0600 Subject: [PATCH 25/26] Raise `NotImplementedError` for `RayVectorArray.volume_cell` named-arrays 2.4.0 added `volume_cell` to `AbstractSpectralPositionalVectorArray`, which `RayVectorArray` subclasses. A ray bundle is a scattered collection of rays rather than a logically-rectangular grid, so the per-voxel volume is undefined and the inherited implementation crashes on scalar-wavelength ray fixtures. Override it to raise `NotImplementedError`, matching the existing `type_matrix`/`test_matrix` precedent and the test's expectation. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/rays/_ray_vectors.py | 8 ++++++++ 1 file changed, 8 insertions(+) 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 From a2003ec5d8bdd6586abfa0ad5e23bf910cb2a512 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Fri, 31 Jul 2026 14:22:25 -0600 Subject: [PATCH 26/26] Address code-review findings on the linear-system PR - backproject_from_weights: re-attach the full target wavelength grid before `photons_absorbed`, so the integrated readout is spread across the correct number of wavelength bins. The detector readout only carries the two band edges, so `num_wavelength` was always 1 and the per-bin spreading (and per-bin quantum efficiency) never happened. - AbstractSystem.backproject: make it a concrete default that raises `NotImplementedError` instead of an `@abc.abstractmethod`, since backprojection is an optional capability (a `SequentialSystem` has no transpose). This stops it from breaking unrelated `AbstractSystem` subclasses that never backproject. - image/backproject and the *_from_weights variants: make the optional arguments keyword-only, so a positional `image(scene, pupil)` can no longer silently bind `pupil` to `axis_wavelength` after the parameter reorder. - _collapse_wavelength: use min/max rather than the first/last edge, so a non-ascending wavelength grid still reports sorted band edges. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ --- optika/sensors/_sensors.py | 4 ++-- optika/systems/_linear.py | 16 ++++++++++++++++ optika/systems/_sequential.py | 2 ++ optika/systems/_systems.py | 11 ++++++++++- 4 files changed, 30 insertions(+), 3 deletions(-) diff --git a/optika/sensors/_sensors.py b/optika/sensors/_sensors.py index aa27105f..9b7000e8 100644 --- a/optika/sensors/_sensors.py +++ b/optika/sensors/_sensors.py @@ -201,8 +201,8 @@ def _collapse_wavelength( return inputs.replace( wavelength=na.stack( arrays=[ - wavelength[{axis_wavelength: +0}], - wavelength[{axis_wavelength: ~0}], + wavelength.min(axis_wavelength), + wavelength.max(axis_wavelength), ], axis=axis_wavelength, ) diff --git a/optika/systems/_linear.py b/optika/systems/_linear.py index 6e2dccba..4380c606 100644 --- a/optika/systems/_linear.py +++ b/optika/systems/_linear.py @@ -272,6 +272,7 @@ 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, @@ -392,6 +393,7 @@ def backproject_from_weights( 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, @@ -469,6 +471,18 @@ def backproject_from_weights( (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( @@ -508,6 +522,7 @@ def backproject_from_weights( 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, @@ -602,6 +617,7 @@ def backproject( | na.AbstractScalar ), coordinates: na.SpectralPositionalVectorArray, + *, axis_wavelength: None | str = None, axis_field: None | tuple[str, str] = None, integrate: bool = True, diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index e548231b..a1e617a6 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -1090,6 +1090,7 @@ def _rayfunction_from_vertices( 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, @@ -1203,6 +1204,7 @@ 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, diff --git a/optika/systems/_systems.py b/optika/systems/_systems.py index fb98289e..19da5dd3 100644 --- a/optika/systems/_systems.py +++ b/optika/systems/_systems.py @@ -27,6 +27,7 @@ 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, @@ -57,11 +58,11 @@ def image( of this method. """ - @abc.abstractmethod 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, @@ -72,6 +73,11 @@ def backproject( 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 @@ -96,3 +102,6 @@ def backproject( Additional keyword arguments used by subclass implementations of this method. """ + raise NotImplementedError( # pragma: nocover + f"{type(self).__name__} does not support backprojection." + )