diff --git a/optika/apertures/_apertures.py b/optika/apertures/_apertures.py index f30acd1..4d45833 100644 --- a/optika/apertures/_apertures.py +++ b/optika/apertures/_apertures.py @@ -600,29 +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: - 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 + center, half = self._bound_center_half() + return center - half @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 + center, half = self._bound_center_half() + return center + half @property def vertices(self) -> None: @@ -699,11 +719,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 +882,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..4cbc4ec 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,67 @@ 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))) + + +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)