Fixed aperture bounding boxes ignoring the internal transformation

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]

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)