From c7186f3dd1a05c53989454b360be09f348f7a790 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 27 Aug 2026 14:52:34 -0600 Subject: [PATCH 1/2] Solve for the stops once when linearizing, not three times Denormalizing a grid does two unrelated things: it solves for where the stops put the field and pupil, and it applies an affine map. The solve costs 0.13 s and depends on nothing but the wavelengths; the map is free. `linearize` denormalizes three grids on the same wavelengths, one per fit, so it solved three times for the same answer, and a fourth for `rayfunction_default`. Split the two halves apart, so that a caller with several grids can solve once: _calc_rayfunction_stops_denormalize the solve _denormalize_grid_from_rays the affine map _denormalize_grid both, as before and have `linearize` solve once, denormalize both of its grids, and pass each fit a grid which is already physical, so their own denormalization is a no-op. before: 4 stop solves, 0.59 s of 3.98 s after: 2 stop solves, 0.29 s of 3.59 s The one which remains is `rayfunction_default`'s, on the default grid rather than the one being linearized. Since the fits can no longer be handed a grid of `None`, `linearize` resolves the defaults itself, and the pupil grid `area_effective` invents when given none becomes `_pupil_vertices_default` so that both can reach it. The two fits keep falling back to `grid_input.pupil` as before: those grids differ because one wants cell vertices and the other wants sample points, and reconciling them is #187's job, not this one's. Every product of `linearize` is unchanged bit for bit: all twenty distortion coefficients, all ten vignetting coefficients, the direction, the scene and sensor coordinates the distortion is fit to, and the illumination the vignetting is fit to. Co-Authored-By: Claude Opus 5 --- optika/systems/_sequential.py | 131 +++++++++++++++++++++++++++++----- 1 file changed, 112 insertions(+), 19 deletions(-) diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index 59381c3..6f736a1 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -781,6 +781,47 @@ def pupil_max(self): """ return self.pupil_boundary.max(self._axis_stops) + @property + def _pupil_vertices_default(self) -> na.Cartesian2dVectorArray: + """ + The pupil grid :meth:`area_effective` uses when given none. + + Its components are cell *vertices*, since the effective area weights + each ray by the area of its pupil cell, whereas + :attr:`grid_input` holds the points that rays are traced at. That is + why the two cannot share a default. + """ + return na.Cartesian2dVectorArray( + x=na.linspace(-1, 1, axis="_pupil_x", num=11), + y=na.linspace(-1, 1, axis="_pupil_y", num=11), + ) + + def _calc_rayfunction_stops_denormalize( + self, + wavelength: na.ScalarLike, + ) -> optika.rays.RayFunctionArray: + """ + Solve for the stops on the sampling that :meth:`_denormalize_grid` + uses, which is the expensive half of denormalizing a grid. + + Separated so that a caller denormalizing several grids on the same + wavelengths can solve once and pass the result to + :meth:`_denormalize_grid_from_rays`, instead of repeating a solve + which depends on nothing else. + + Parameters + ---------- + wavelength + The wavelengths to solve at. + """ + return self._calc_rayfunction_stops( + wavelength_input=wavelength, + axis_pupil_stop=self._axis_pupil_stop, + axis_field_stop=self._axis_field_stop, + samples_pupil_stop=21, + samples_field_stop=21, + ) + def _denormalize_grid( self, grid: optika.vectors.ObjectVectorArray, @@ -791,17 +832,41 @@ def _denormalize_grid( if (not normalized_field) and (not normalized_pupil): return grid + return self._denormalize_grid_from_rays( + grid=grid, + rayfunction_stops=self._calc_rayfunction_stops_denormalize(grid.wavelength), + normalized_field=normalized_field, + normalized_pupil=normalized_pupil, + ) + + def _denormalize_grid_from_rays( + self, + grid: optika.vectors.ObjectVectorArray, + rayfunction_stops: optika.rays.RayFunctionArray, + normalized_field: bool = True, + normalized_pupil: bool = True, + ) -> optika.vectors.ObjectVectorArray: + """ + Map a normalized grid onto the field and pupil that the given stop + rayfunction describes. + + Parameters + ---------- + grid + The grid to denormalize. + rayfunction_stops + The result of :meth:`_calc_rayfunction_stops_denormalize` on the + wavelengths of `grid`. + normalized_field + A boolean flag indicating whether the field of `grid` is given in + normalized or physical units. + normalized_pupil + A boolean flag indicating whether the pupil of `grid` is given in + normalized or physical units. + """ axis_field = self._axis_field_stop axis_pupil = self._axis_pupil_stop - rayfunction_stops = self._calc_rayfunction_stops( - wavelength_input=grid.wavelength, - axis_pupil_stop=axis_pupil, - axis_field_stop=axis_field, - samples_pupil_stop=21, - samples_field_stop=21, - ) - result = grid.copy_shallow() object_is_at_infinity = self.object_is_at_infinity @@ -1526,10 +1591,7 @@ def area_effective( if field is None: field = self.grid_input.field if pupil is None: - pupil = na.Cartesian2dVectorArray( - x=na.linspace(-1, 1, axis="_pupil_x", num=11), - y=na.linspace(-1, 1, axis="_pupil_y", num=11), - ) + pupil = self._pupil_vertices_default grid = optika.vectors.ObjectVectorArray( wavelength=wavelength, @@ -1665,16 +1727,47 @@ def linearize( # `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. + # With no pupil given the two fall back to different grids, since + # `grid_input.pupil` holds points and cannot serve as vertices. if pupil is not None: pupil_centers = pupil.cell_centers(axis=tuple(na.shape(pupil))) else: - pupil_centers = None + pupil = self._pupil_vertices_default + pupil_centers = self.grid_input.pupil + + if wavelength is None: + wavelength = self.grid_input.wavelength + if field is None: + field = self.grid_input.field + + # Each of the three fits below denormalizes its own grid, and the + # expensive part of that is solving for the stops, which depends on + # nothing but the wavelengths. Solve once here and hand each of them + # a grid which is already physical. + stops = self._calc_rayfunction_stops_denormalize(wavelength) + + def denormalize( + pupil: na.AbstractCartesian2dVectorArray, + ) -> optika.vectors.ObjectVectorArray: + return self._denormalize_grid_from_rays( + grid=optika.vectors.ObjectVectorArray( + wavelength=wavelength, + field=field, + pupil=pupil, + ), + rayfunction_stops=stops, + normalized_field=normalized_field, + normalized_pupil=normalized_pupil, + ) + + grid_area = denormalize(pupil) + grid_fit = denormalize(pupil_centers) kwargs = dict( wavelength=wavelength, - field=field, - normalized_field=normalized_field, - normalized_pupil=normalized_pupil, + field=grid_area.field, + normalized_field=False, + normalized_pupil=False, ) # the cosine of the refracted angle at which light strikes the sensor, # computed the same way as @@ -1693,11 +1786,11 @@ def linearize( ) return LinearSystem( - area_effective=self.area_effective(pupil=pupil, **kwargs), - distortion=self.distortion(pupil=pupil_centers, degree=degree, **kwargs), + area_effective=self.area_effective(pupil=grid_area.pupil, **kwargs), + distortion=self.distortion(pupil=grid_fit.pupil, degree=degree, **kwargs), sensor=self.sensor, direction=direction, - vignetting=self.vignetting(pupil=pupil_centers, degree=degree, **kwargs), + vignetting=self.vignetting(pupil=grid_fit.pupil, degree=degree, **kwargs), ) def _rayfunction_and_axes( From 0e0f0f11835fc5769f9f2d5d31a722e447438830 Mon Sep 17 00:00:00 2001 From: Roy Smart Date: Thu, 27 Aug 2026 16:06:03 -0600 Subject: [PATCH 2/2] Drop the wrapper by fixing the solver's defaults `_calc_rayfunction_stops_denormalize` existed only to fill in four arguments which every caller passes identically. It is not needed: the defaults on `_calc_rayfunction_stops` were simply wrong. Its axes had no defaults at all even though both callers pass the class constants, and its sample counts defaulted to 101, which nothing has ever used. Give it the values its callers actually want, and the wrapper and the boilerplate in `rayfunction_stops` both go away. Net 27 lines fewer, and the products of `linearize` remain identical to main across all 37 arrays. --- optika/systems/_sequential.py | 55 +++++++++-------------------------- 1 file changed, 14 insertions(+), 41 deletions(-) diff --git a/optika/systems/_sequential.py b/optika/systems/_sequential.py index 6f736a1..f4683fd 100644 --- a/optika/systems/_sequential.py +++ b/optika/systems/_sequential.py @@ -636,11 +636,16 @@ def jacobian(x, _function=function, _dx=dx): def _calc_rayfunction_stops( self, wavelength_input: na.ScalarLike, - axis_pupil_stop: str, - axis_field_stop: str, - samples_pupil_stop: int = 101, - samples_field_stop: int = 101, + axis_pupil_stop: None | str = None, + axis_field_stop: None | str = None, + samples_pupil_stop: int = 21, + samples_field_stop: int = 21, ) -> optika.rays.RayFunctionArray: + if axis_pupil_stop is None: + axis_pupil_stop = self._axis_pupil_stop + if axis_field_stop is None: + axis_field_stop = self._axis_field_stop + surfaces = self.surfaces_all index_pupil_stop = self.index_pupil_stop @@ -699,13 +704,7 @@ def rayfunction_stops(self) -> optika.rays.RayFunctionArray: which is designed to exactly strike the borders of both the field stop and the pupil stop. """ - return self._calc_rayfunction_stops( - wavelength_input=self.grid_input.wavelength, - axis_pupil_stop=self._axis_pupil_stop, - axis_field_stop=self._axis_field_stop, - samples_pupil_stop=21, - samples_field_stop=21, - ) + return self._calc_rayfunction_stops(self.grid_input.wavelength) @property def _axis_stops(self) -> tuple[str, str]: @@ -796,32 +795,6 @@ def _pupil_vertices_default(self) -> na.Cartesian2dVectorArray: y=na.linspace(-1, 1, axis="_pupil_y", num=11), ) - def _calc_rayfunction_stops_denormalize( - self, - wavelength: na.ScalarLike, - ) -> optika.rays.RayFunctionArray: - """ - Solve for the stops on the sampling that :meth:`_denormalize_grid` - uses, which is the expensive half of denormalizing a grid. - - Separated so that a caller denormalizing several grids on the same - wavelengths can solve once and pass the result to - :meth:`_denormalize_grid_from_rays`, instead of repeating a solve - which depends on nothing else. - - Parameters - ---------- - wavelength - The wavelengths to solve at. - """ - return self._calc_rayfunction_stops( - wavelength_input=wavelength, - axis_pupil_stop=self._axis_pupil_stop, - axis_field_stop=self._axis_field_stop, - samples_pupil_stop=21, - samples_field_stop=21, - ) - def _denormalize_grid( self, grid: optika.vectors.ObjectVectorArray, @@ -834,7 +807,7 @@ def _denormalize_grid( return self._denormalize_grid_from_rays( grid=grid, - rayfunction_stops=self._calc_rayfunction_stops_denormalize(grid.wavelength), + rayfunction_stops=self._calc_rayfunction_stops(grid.wavelength), normalized_field=normalized_field, normalized_pupil=normalized_pupil, ) @@ -855,8 +828,8 @@ def _denormalize_grid_from_rays( grid The grid to denormalize. rayfunction_stops - The result of :meth:`_calc_rayfunction_stops_denormalize` on the - wavelengths of `grid`. + The result of :meth:`_calc_rayfunction_stops` on the wavelengths + of `grid`. normalized_field A boolean flag indicating whether the field of `grid` is given in normalized or physical units. @@ -1744,7 +1717,7 @@ def linearize( # expensive part of that is solving for the stops, which depends on # nothing but the wavelengths. Solve once here and hand each of them # a grid which is already physical. - stops = self._calc_rayfunction_stops_denormalize(wavelength) + stops = self._calc_rayfunction_stops(wavelength) def denormalize( pupil: na.AbstractCartesian2dVectorArray,