From 1cefb0ded14d30ca876abfea568ad7404f0f72a3 Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Thu, 11 Jun 2026 10:25:32 -0600 Subject: [PATCH 1/2] Fix aperture bounding boxes ignoring the internal transformation AbstractPolygonalAperture.bound_lower/bound_upper were computed from the raw vertices, while __call__ and wire() apply the aperture's internal transformation. For apertures with a nonzero internal transformation (e.g. a decentered grating aperture), any consumer of the bounding box (stratified wavefield sampling, sensor pixel grids, pupil normalization) saw the wrong region of the surface. EllipticalAperture had the related defect of transforming only the bounding-box corner points, which is incorrect under rotation; its bounds are now derived from the wire, matching CircularSectorAperture. RectangularAperture.__call__ expressed the rectangle in terms of bound_lower/bound_upper while also inverse-transforming the position, which double-counted the internal transformation once the bounds included it; it now tests against half_width directly. The generic bound tests now assert that the bounds enclose the wire and that the wire centroid is accepted by active, non-inverted apertures (nominal values, with a small tolerance for interpolation round-off, since uncertain parameters may be redrawn between evaluations). Co-Authored-By: Claude Fable 5 --- optika/apertures/_apertures.py | 46 ++++++++------------- optika/apertures/_apertures_test.py | 64 +++++++++++++++++++++++++++-- 2 files changed, 78 insertions(+), 32 deletions(-) diff --git a/optika/apertures/_apertures.py b/optika/apertures/_apertures.py index f30acd1..45f4f84 100644 --- a/optika/apertures/_apertures.py +++ b/optika/apertures/_apertures.py @@ -602,27 +602,11 @@ def __call__( @property def bound_lower(self) -> na.Cartesian3dVectorArray: - unit = na.unit(self.radius) - result = na.Cartesian3dVectorArray() - if unit is not None: - result = result * unit - if self.transformation is not None: - result = self.transformation(result) - result.x = result.x - self.radius.x - result.y = result.y - self.radius.y - return result + return self.wire().min() @property def bound_upper(self) -> na.Cartesian3dVectorArray: - unit = na.unit(self.radius) - result = na.Cartesian3dVectorArray() - if unit is not None: - result = result * unit - if self.transformation is not None: - result = self.transformation(result) - result.x = result.x + self.radius.x - result.y = result.y + self.radius.y - return result + return self.wire().max() @property def vertices(self) -> None: @@ -699,11 +683,17 @@ def __call__( @property def bound_lower(self) -> na.AbstractCartesian3dVectorArray: - return self.vertices.min(axis="vertex") + vertices = self.vertices + if self.transformation is not None: + vertices = self.transformation(vertices) + return vertices.min(axis="vertex") @property def bound_upper(self) -> na.AbstractCartesian3dVectorArray: - return self.vertices.max(axis="vertex") + vertices = self.vertices + if self.transformation is not None: + vertices = self.transformation(vertices) + return vertices.max(axis="vertex") def wire(self, num: None | int = None) -> na.Cartesian3dVectorArray: if num is None: @@ -856,24 +846,24 @@ def __call__( self, position: na.AbstractCartesian3dVectorArray, ) -> na.AbstractScalar: - bound_lower = self.bound_lower - bound_upper = self.bound_upper + half_width = na.asanyarray( + self.half_width, + like=na.Cartesian2dVectorArray(), + ) active = self.active inverted = self.inverted if self.transformation is not None: position = self.transformation.inverse(position) + position = position.xy - shape = na.shape_broadcasted( - bound_lower, bound_upper, active, inverted, position - ) + shape = na.shape_broadcasted(half_width, active, inverted, position) - bound_lower = na.broadcast_to(bound_lower, shape) - bound_upper = na.broadcast_to(bound_upper, shape) + half_width = na.broadcast_to(half_width, shape) active = na.broadcast_to(active, shape) inverted = na.broadcast_to(inverted, shape) position = na.broadcast_to(position, shape) - mask = (bound_lower <= position) & (position <= bound_upper) + mask = (-half_width <= position) & (position <= half_width) mask = mask.x & mask.y mask[inverted] = ~mask[inverted] diff --git a/optika/apertures/_apertures_test.py b/optika/apertures/_apertures_test.py index c6ae0a8..a02ea51 100644 --- a/optika/apertures/_apertures_test.py +++ b/optika/apertures/_apertures_test.py @@ -25,6 +25,12 @@ ] +def _nominal(x: na.AbstractScalar) -> na.AbstractScalar: + if isinstance(x, na.AbstractUncertainScalarArray): + return x.nominal + return x + + class AbstractTestAbstractAperture( test_mixins.AbstractTestDxfWritable, test_mixins.AbstractTestPrintable, @@ -57,6 +63,17 @@ def test__call__(self, a: optika.apertures.AbstractAperture): assert np.all(result[~na.as_named_array(a.active)]) + # the centroid of the wire must be inside an active, non-inverted + # aperture (every aperture here is star-shaped about its centroid) + if (a.active is True) and (a.inverted is False): + centroid = a.wire().mean("wire") + centroid = na.Cartesian3dVectorArray( + x=_nominal(centroid.x), + y=_nominal(centroid.y), + z=_nominal(centroid.z), + ) + assert np.all(_nominal(na.as_named_array(a(centroid)))) + @pytest.mark.parametrize("rays", optika.rays._tests.test_ray_vectors.rays) def test_clip_rays( self, @@ -75,12 +92,28 @@ def test_clip_rays( assert np.mean(result.unvignetted) <= np.mean(rays.unvignetted) def test_bound_lower(self, a: optika.apertures.AbstractAperture): - assert isinstance(a.bound_lower, na.AbstractCartesian3dVectorArray) + result = a.bound_lower + assert isinstance(result, na.AbstractCartesian3dVectorArray) + # compare nominal values only, since uncertain parameters may be + # redrawn between independent evaluations of the geometry + wire = a.wire() + for component in ("x", "y"): + bound = _nominal(getattr(result, component)) + edge = _nominal(getattr(wire, component).min("wire")) + tolerance = 1e-9 * (np.abs(bound) + np.abs(edge)) + assert np.all(bound <= edge + tolerance) def test_bound_upper(self, a: optika.apertures.AbstractAperture): - assert isinstance(a.bound_upper, na.AbstractCartesian3dVectorArray) - assert np.all(a.bound_upper.x != a.bound_lower.x) - assert np.all(a.bound_upper.y != a.bound_lower.y) + result = a.bound_upper + assert isinstance(result, na.AbstractCartesian3dVectorArray) + assert np.all(result.x != a.bound_lower.x) + assert np.all(result.y != a.bound_lower.y) + wire = a.wire() + for component in ("x", "y"): + bound = _nominal(getattr(result, component)) + edge = _nominal(getattr(wire, component).max("wire")) + tolerance = 1e-9 * (np.abs(bound) + np.abs(edge)) + assert np.all(bound >= edge - tolerance) def test_vertices(self, a: optika.apertures.AbstractAperture): if a.vertices is not None: @@ -386,3 +419,26 @@ class TestIsoscelesTrapezoidalAperture( AbstractTestAbstractIsoscelesTrapezoidalAperture, ): pass + + +def test_rectangular_aperture_decentered(): + """ + A rectangular aperture decentered by its internal transformation must + accept points inside the translated rectangle and reject points inside + the untranslated one. + + Regression test: ``RectangularAperture.__call__`` used to express the + rectangle in terms of :attr:`bound_lower`/:attr:`bound_upper` while + also inverse-transforming the position, applying the internal + transformation twice. + """ + half_width = 5 * u.mm + decenter = 12 * u.mm + a = optika.apertures.RectangularAperture( + half_width=half_width, + transformation=na.transformations.Cartesian3dTranslation(x=decenter), + ) + point_inside = na.Cartesian3dVectorArray(12, 0, 0) * u.mm + point_outside = na.Cartesian3dVectorArray(0, 0, 0) * u.mm + assert np.all(na.as_named_array(a(point_inside))) + assert not np.any(na.as_named_array(a(point_outside))) From 947977119b7aa5a1526c746c89b1834fac2ab601 Mon Sep 17 00:00:00 2001 From: jacobdparker Date: Mon, 15 Jun 2026 11:42:38 -0600 Subject: [PATCH 2/2] Compute elliptical aperture bounds analytically Address review feedback: deriving EllipticalAperture.bound_lower/upper from wire().min()/max() sampled the boundary at a finite number of points, so the true extent between samples was underestimated whenever the internal transformation rotated the ellipse off the coordinate axes. Replace the wire sampling with a closed form. A boundary point is p(t) = c + a*e_a*cos(t) + b*e_b*sin(t), where c is the transformed center and e_a, e_b are the transformed local axes, so the extent along any world component is sqrt((a*e_a)^2 + (b*e_b)^2). This is exact under rotation and reduces to the axis-aligned (radius.x, radius.y) box when untransformed. CircularAperture bounds are already analytic and rotation-invariant, so they are unchanged. Co-Authored-By: Claude Opus 4.8 --- optika/apertures/_apertures.py | 40 ++++++++++++++++++++++++++-- optika/apertures/_apertures_test.py | 41 +++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+), 2 deletions(-) diff --git a/optika/apertures/_apertures.py b/optika/apertures/_apertures.py index 45f4f84..4d45833 100644 --- a/optika/apertures/_apertures.py +++ b/optika/apertures/_apertures.py @@ -600,13 +600,49 @@ def __call__( return mask + def _bound_center_half( + self, + ) -> tuple[na.Cartesian3dVectorArray, na.Cartesian3dVectorArray]: + """ + The center and per-component half-extent of the axis-aligned bounding + box, computed analytically so the bound is exact even when + :attr:`transformation` rotates the ellipse. + + A point on the ellipse boundary is + :math:`p(t) = c + a\\,\\hat{e}_a \\cos t + b\\,\\hat{e}_b \\sin t`, + where :math:`c` is the center and :math:`\\hat{e}_a, \\hat{e}_b` are the + images of the local axes under the transformation. The extent along + any world component is then + :math:`\\sqrt{(a\\,\\hat{e}_a)^2 + (b\\,\\hat{e}_b)^2}`, since + :math:`\\max_t (A \\cos t + B \\sin t) = \\sqrt{A^2 + B^2}`. + """ + radius = self.radius + zero = na.Cartesian3dVectorArray() * radius.x + axis_a = na.Cartesian3dVectorArray(x=radius.x, y=0 * radius.x, z=0 * radius.x) + axis_b = na.Cartesian3dVectorArray(x=0 * radius.y, y=radius.y, z=0 * radius.y) + + center = zero + if self.transformation is not None: + center = self.transformation(zero) + axis_a = self.transformation(axis_a) - center + axis_b = self.transformation(axis_b) - center + + half = na.Cartesian3dVectorArray( + x=np.sqrt(np.square(axis_a.x) + np.square(axis_b.x)), + y=np.sqrt(np.square(axis_a.y) + np.square(axis_b.y)), + z=np.sqrt(np.square(axis_a.z) + np.square(axis_b.z)), + ) + return center, half + @property def bound_lower(self) -> na.Cartesian3dVectorArray: - return self.wire().min() + center, half = self._bound_center_half() + return center - half @property def bound_upper(self) -> na.Cartesian3dVectorArray: - return self.wire().max() + center, half = self._bound_center_half() + return center + half @property def vertices(self) -> None: diff --git a/optika/apertures/_apertures_test.py b/optika/apertures/_apertures_test.py index a02ea51..4cbc4ec 100644 --- a/optika/apertures/_apertures_test.py +++ b/optika/apertures/_apertures_test.py @@ -442,3 +442,44 @@ def test_rectangular_aperture_decentered(): point_outside = na.Cartesian3dVectorArray(0, 0, 0) * u.mm assert np.all(na.as_named_array(a(point_inside))) assert not np.any(na.as_named_array(a(point_outside))) + + +def test_elliptical_aperture_bounds_analytic(): + """ + The bounding box of a rotated, decentered elliptical aperture is computed + analytically and must match the true extent of the ellipse. + + Regression test: the bounds used to be sampled from a coarse + :meth:`wire`, which underestimates the extent between samples whenever the + transformation rotates the ellipse off the coordinate axes. + """ + a = optika.apertures.EllipticalAperture( + radius=na.Cartesian2dVectorArray(100, 50) * u.mm, + transformation=( + na.transformations.Cartesian3dRotationZ(30 * u.deg) + @ na.transformations.Cartesian3dTranslation(x=5 * u.mm, y=-2 * u.mm) + ), + ) + + bound_lower = a.bound_lower + bound_upper = a.bound_upper + + # a densely-sampled wire approaches the true extent from the inside + wire = a.wire(num=100001) + for component in ("x", "y"): + lower = getattr(bound_lower, component) + upper = getattr(bound_upper, component) + edge_lower = getattr(wire, component).min("wire") + edge_upper = getattr(wire, component).max("wire") + scale = np.abs(upper - lower) + + # the analytic bound encloses the dense wire ... + assert np.all(lower <= edge_lower + 1e-9 * scale) + assert np.all(upper >= edge_upper - 1e-9 * scale) + # ... and is tight: it does not overshoot it + assert np.all(np.abs(lower - edge_lower) < 1e-6 * scale) + assert np.all(np.abs(upper - edge_upper) < 1e-6 * scale) + + # an in-plane rotation leaves the aperture flat in z + assert np.all(bound_lower.z == 0 * u.mm) + assert np.all(bound_upper.z == 0 * u.mm)