From b1f88f4a7424a593915dabd9b00f2ed0de256a74 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 2 Sep 2026 17:57:49 -0400 Subject: [PATCH 1/3] feat: Shape.contains and Shape.boundary for source-plane shapes (image-source-mappings phase 2a) Point containment and closed boundary polygons for Circle, Triangle, Polygon and Square, documenting the (y, x) axis convention the triangle solvers already use. Fixes the reflected containment test in Triangle.triangle_contains_mask, which unpacked its vertices as (y, x) while testing a centroid built as (x, y). Square.contains / Square.boundary sort their bounds so they work for either top/bottom ordering; mask and area keep their documented assumption. Part of PyAutoLabs/PyAutoLens#719 (image-source-mappings epic, phase 2a). Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_011EVnQ4sEYrRM7GCvFyveUA --- autoarray/structures/triangles/shape.py | 396 +++++++++++++++++- .../structures/triangles/test_shape.py | 184 ++++++++ 2 files changed, 559 insertions(+), 21 deletions(-) create mode 100644 test_autoarray/structures/triangles/test_shape.py diff --git a/autoarray/structures/triangles/shape.py b/autoarray/structures/triangles/shape.py index cba3056bf..c58b0adcf 100644 --- a/autoarray/structures/triangles/shape.py +++ b/autoarray/structures/triangles/shape.py @@ -8,6 +8,26 @@ class Shape(ABC): """ A shape in the source plane for which we identify corresponding image plane pixels using up-sampling. + + Coordinate convention + --------------------- + Coordinates are indexed here exactly as they are in the triangle arrays: + element ``0`` of a coordinate pair is the attribute these classes call + ``x`` and element ``1`` is the attribute they call ``y``. For a triangle + array of shape ``(n, 3, 2)`` this is ``triangles[..., 0]`` and + ``triangles[..., 1]``, and the ``(N, 2)`` array of points passed to + `contains` is ordered the same way, as is the polygon returned by + `boundary`. + + The legacy attribute names do **not** correspond to the physical axes. + Every grid PyAuto produces is ordered ``(y, x)``, and ``PointSolver.solve`` + builds its shape as ``Point(*source_plane_coordinate)`` from such a + ``(y, x)`` tuple, so the attribute named `x` in fact holds the ``y`` + (first) coordinate and the attribute named `y` holds the ``x`` (second) + coordinate. The names are kept for backwards compatibility; what matters + is that points, vertices and triangle arrays all use element ``0`` for the + first axis, so containment is computed consistently and a `boundary` can be + plotted directly against an arcsec ``(y, x)`` grid. """ @property @@ -32,6 +52,92 @@ def mask(self, triangles: np.ndarray) -> np.ndarray: A boolean array indicating which triangles contain the shape. """ + @abstractmethod + def contains(self, points: np.ndarray) -> np.ndarray: + """ + Determine which points lie inside the shape. + + Parameters + ---------- + points + An array of coordinates of shape ``(N, 2)``, ordered the same way + as the triangle vertices (see the class docstring). + + Returns + ------- + A boolean array of shape ``(N,)`` indicating which points lie inside + the shape. + """ + + @abstractmethod + def boundary(self, n: int = 100) -> np.ndarray: + """ + The boundary of the shape, as a closed polygon. + + Parameters + ---------- + n + The number of samples used for curved boundaries (e.g. a circle). + Ignored by shapes whose boundary is already polygonal. + + Returns + ------- + An array of shape ``(M, 2)`` whose first row is repeated as its last + row, ordered the same way as the triangle vertices. + """ + + +def _barycentric_contains( + a_0, + a_1, + b_0, + b_1, + c_0, + c_1, + coordinate_0, + coordinate_1, +) -> np.ndarray: + """ + Determine which coordinates lie inside the triangle with vertices + ``(a_0, a_1)``, ``(b_0, b_1)`` and ``(c_0, c_1)``, using barycentric + coordinates. + + Vertices and coordinates are indexed with the same axis order (element 0 + against element 0) — see the `Shape` docstring. Either the vertices or the + coordinates may be arrays, so this single implementation serves both + "which triangles contain this shape" and "which points does this shape + contain". + + Parameters + ---------- + a_0, a_1, b_0, b_1, c_0, c_1 + The components of the three vertices of the triangle. + coordinate_0, coordinate_1 + The components of the coordinates being tested. + + Returns + ------- + A boolean array indicating which coordinates lie inside the triangle. + """ + denominator = (b_1 - c_1) * (a_0 - c_0) + (c_0 - b_0) * (a_1 - c_1) + + alpha = ( + (b_1 - c_1) * (coordinate_0 - c_0) + (c_0 - b_0) * (coordinate_1 - c_1) + ) / denominator + beta = ( + (c_1 - a_1) * (coordinate_0 - c_0) + (a_0 - c_0) * (coordinate_1 - c_1) + ) / denominator + gamma = 1 - alpha - beta + + return ( + (0 <= alpha) + & (alpha <= 1) + & (0 <= beta) + & (beta <= 1) + & (0 <= gamma) + & (gamma <= 1) + ) + class Point(Shape): def __init__(self, x: float, y: float): @@ -68,17 +174,49 @@ def mask(self, triangles: np.ndarray) -> np.ndarray: ------- A boolean array indicating which triangles contain the point. """ - y1, x1 = triangles[:, 0, 1], triangles[:, 0, 0] - y2, x2 = triangles[:, 1, 1], triangles[:, 1, 0] - y3, x3 = triangles[:, 2, 1], triangles[:, 2, 0] + return _barycentric_contains( + triangles[:, 0, 0], + triangles[:, 0, 1], + triangles[:, 1, 0], + triangles[:, 1, 1], + triangles[:, 2, 0], + triangles[:, 2, 1], + self.x, + self.y, + ) - denominator = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3) + def contains(self, points: np.ndarray) -> np.ndarray: + """ + A point has no area, so no coordinate lies inside it. - a = ((y2 - y3) * (self.x - x3) + (x3 - x2) * (self.y - y3)) / denominator - b = ((y3 - y1) * (self.x - x3) + (x1 - x3) * (self.y - y3)) / denominator - c = 1 - a - b + Raises + ------ + NotImplementedError + Always. Use `Circle` (or another finite `Shape`) for a source-plane + region with an interior, or `PointSolver` to find the image-plane + positions to which a point source traces. + """ + raise NotImplementedError( + "A Point has zero area so it contains no coordinates. Use a Circle " + "(or another finite Shape) for a source-plane region with an " + "interior, or PointSolver to find the image plane positions of a " + "point source." + ) + + def boundary(self, n: int = 100) -> np.ndarray: + """ + The boundary of a point is the point itself, a single row. - return (0 <= a) & (a <= 1) & (0 <= b) & (b <= 1) & (0 <= c) & (c <= 1) + Parameters + ---------- + n + Unused; a point is not sampled. + + Returns + ------- + An array of shape ``(1, 2)`` containing the coordinates of the point. + """ + return np.array([[self.x, self.y]]) def tree_flatten(self): """ @@ -161,6 +299,53 @@ def mask(self, triangles: np.ndarray) -> np.ndarray: return (distance_squared <= radius_2) | super().mask(triangles) + def contains(self, points: np.ndarray) -> np.ndarray: + """ + Determine which points lie inside the circle. + + Parameters + ---------- + points + An array of coordinates of shape ``(N, 2)``, ordered the same way + as the triangle vertices (see the `Shape` docstring). + + Returns + ------- + A boolean array of shape ``(N,)``; points exactly on the circle count + as inside. + """ + points = np.asarray(points) + + delta_0 = points[:, 0] - self.x + delta_1 = points[:, 1] - self.y + + return delta_0 * delta_0 + delta_1 * delta_1 <= self.radius * self.radius + + def boundary(self, n: int = 100) -> np.ndarray: + """ + The circle sampled as a closed polygon of ``n`` points. + + Parameters + ---------- + n + The number of samples around the circle. + + Returns + ------- + An array of shape ``(n + 1, 2)`` whose last row repeats its first. + """ + angles = np.linspace(0.0, 2.0 * np.pi, n, endpoint=False) + + boundary = np.stack( + ( + self.x + self.radius * np.cos(angles), + self.y + self.radius * np.sin(angles), + ), + axis=-1, + ) + + return np.append(boundary, boundary[:1], axis=0) + def tree_flatten(self): """ Flatten this model as a PyTree. @@ -227,23 +412,80 @@ def mask(self, triangles: np.ndarray) -> np.ndarray: return self.triangle_contains_mask(triangles) | super().mask(triangles) def triangle_contains_mask(self, triangles: np.ndarray) -> np.ndarray: - y1, x1 = self.a - y2, x2 = self.b - y3, x3 = self.c + """ + Determine which triangles have their centroid inside this triangle. - denominator = (y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3) + The vertices of this triangle and the centroids being tested are + indexed with the same axis order (element 0 against element 0). This + previously unpacked the vertices as ``y, x`` while testing a centroid + built as ``x, y``, which tested the reflected triangle and so gave the + wrong answer for any triangle that is not symmetric about + ``element 0 == element 1``. - centroid_x, centroid_y = centroid(triangles) + Parameters + ---------- + triangles + The vertices of the triangles. + + Returns + ------- + A boolean array indicating which triangle centroids lie inside this + triangle. + """ + centroid_0, centroid_1 = centroid(triangles) + + return _barycentric_contains( + self.a[0], + self.a[1], + self.b[0], + self.b[1], + self.c[0], + self.c[1], + centroid_0, + centroid_1, + ) + + def contains(self, points: np.ndarray) -> np.ndarray: + """ + Determine which points lie inside the triangle. + + Parameters + ---------- + points + An array of coordinates of shape ``(N, 2)``, ordered the same way + as the vertices (see the `Shape` docstring). + + Returns + ------- + A boolean array of shape ``(N,)``; points on an edge count as inside. + """ + points = np.asarray(points) + + return _barycentric_contains( + self.a[0], + self.a[1], + self.b[0], + self.b[1], + self.c[0], + self.c[1], + points[:, 0], + points[:, 1], + ) - a = ( - (y2 - y3) * (centroid_x - x3) + (x3 - x2) * (centroid_y - y3) - ) / denominator - b = ( - (y3 - y1) * (centroid_x - x3) + (x1 - x3) * (centroid_y - y3) - ) / denominator - c = 1 - a - b + def boundary(self, n: int = 100) -> np.ndarray: + """ + The three vertices of the triangle, closed by repeating the first. - return (0 <= a) & (a <= 1) & (0 <= b) & (b <= 1) & (0 <= c) & (c <= 1) + Parameters + ---------- + n + Unused; a triangle's boundary is exactly polygonal. + + Returns + ------- + An array of shape ``(4, 2)``. + """ + return np.array([self.a, self.b, self.c, self.a]) @property def area(self) -> float: @@ -331,6 +573,47 @@ def mask(self, triangles: np.ndarray) -> np.ndarray: axis=0, ) | super().mask(triangles) + def contains(self, points: np.ndarray) -> np.ndarray: + """ + Determine which points lie inside the polygon. + + The polygon is decomposed into a fan of triangles about its first + vertex, so a point is inside when it is inside any of them. This is + exact for **convex** polygons only; for a concave polygon the fan + covers its convex hull. + + Parameters + ---------- + points + An array of coordinates of shape ``(N, 2)``, ordered the same way + as the vertices (see the `Shape` docstring). + + Returns + ------- + A boolean array of shape ``(N,)``; points on an edge count as inside. + """ + points = np.asarray(points) + + return np.any( + [triangle.contains(points) for triangle in self.triangles], + axis=0, + ) + + def boundary(self, n: int = 100) -> np.ndarray: + """ + The vertices of the polygon, closed by repeating the first. + + Parameters + ---------- + n + Unused; a polygon's boundary is exactly polygonal. + + Returns + ------- + An array of shape ``(len(vertices) + 1, 2)``. + """ + return np.array(list(self.vertices) + [self.vertices[0]]) + class Square(Point): def __init__(self, top, bottom, left, right): @@ -386,3 +669,74 @@ def mask(self, triangles: np.ndarray) -> np.ndarray: & (self.bottom >= centroid_y) & (centroid_y >= self.top) ) | super().mask(triangles) + + @property + def _bounds(self): + """ + The bounds of the square as ``(low_0, high_0, low_1, high_1)``. + + `mask` and `area` assume ``left < right`` and ``top < bottom`` + numerically (coordinates measured from the top left corner of the + image, so `top` is the smaller number). `contains` and `boundary` do + not: they sort the pair, so they behave the same way when a square is + built from arcsec ``(y, x)`` coordinates, where the top of the image is + the *larger* first coordinate. + """ + return ( + min(self.left, self.right), + max(self.left, self.right), + min(self.top, self.bottom), + max(self.top, self.bottom), + ) + + def contains(self, points: np.ndarray) -> np.ndarray: + """ + Determine which points lie inside the square. + + Parameters + ---------- + points + An array of coordinates of shape ``(N, 2)``, ordered the same way + as the triangle vertices (see the `Shape` docstring), so element 0 + is bounded by `left` / `right` and element 1 by `top` / `bottom` — + the same pairing `mask` uses. + + Returns + ------- + A boolean array of shape ``(N,)``; points on an edge count as inside. + """ + points = np.asarray(points) + + low_0, high_0, low_1, high_1 = self._bounds + + return ( + (low_0 <= points[:, 0]) + & (points[:, 0] <= high_0) + & (low_1 <= points[:, 1]) + & (points[:, 1] <= high_1) + ) + + def boundary(self, n: int = 100) -> np.ndarray: + """ + The four corners of the square, closed by repeating the first. + + Parameters + ---------- + n + Unused; a square's boundary is exactly polygonal. + + Returns + ------- + An array of shape ``(5, 2)``. + """ + low_0, high_0, low_1, high_1 = self._bounds + + return np.array( + [ + [low_0, low_1], + [low_0, high_1], + [high_0, high_1], + [high_0, low_1], + [low_0, low_1], + ] + ) diff --git a/test_autoarray/structures/triangles/test_shape.py b/test_autoarray/structures/triangles/test_shape.py new file mode 100644 index 000000000..e64368cca --- /dev/null +++ b/test_autoarray/structures/triangles/test_shape.py @@ -0,0 +1,184 @@ +import numpy as np +import pytest + +from autoarray.structures.triangles.shape import ( + Circle, + Point, + Polygon, + Square, + Triangle, +) + + +def triangle_array_with_centroid(coordinate, half_width=1.0e-3): + """ + A single, non-degenerate triangle whose centroid is `coordinate`. + + Used to compare a `contains` test on a coordinate against the `mask` test + on a triangle sitting at that coordinate. + """ + coordinate_0, coordinate_1 = coordinate + + return np.array( + [ + [ + [coordinate_0, coordinate_1 + 2.0 * half_width], + [coordinate_0 - np.sqrt(3.0) * half_width, coordinate_1 - half_width], + [coordinate_0 + np.sqrt(3.0) * half_width, coordinate_1 - half_width], + ] + ] + ) + + +""" +An asymmetric triangle, written as the (y, x) coordinates every PyAuto grid +uses: it spans 3 in y and 1 in x, so swapping a coordinate's two elements +changes the answer. Every convention test below relies on that asymmetry. +""" +ASYMMETRIC_TRIANGLE_YX = ((0.0, 0.0), (3.0, 0.0), (0.0, 1.0)) +POINT_INSIDE_YX = (2.0, 0.2) + + +def test_contains_uses_the_same_axis_order_as_the_triangle_array(): + """ + `contains` compares element 0 of a point with `Shape.x` and element 1 with + `Shape.y` — the same pairing `mask` uses when it reads `triangles[..., 0]` + and `triangles[..., 1]`. Since every grid PyAuto produces is ordered + `(y, x)`, and `PointSolver.solve` builds `Point(*source_plane_coordinate)` + from such a tuple, the attribute named `x` holds the y coordinate. + """ + triangles = np.array(ASYMMETRIC_TRIANGLE_YX)[None] + + point = POINT_INSIDE_YX + swapped = point[::-1] + + assert Point(*point).mask(triangles) == np.array([True]) + assert Point(*swapped).mask(triangles) == np.array([False]) + + assert Circle(*point, radius=1.0e-8).contains(np.array([point])) == np.array([True]) + assert Circle(*point, radius=1.0e-8).contains(np.array([swapped])) == np.array( + [False] + ) + + triangle = Triangle(*ASYMMETRIC_TRIANGLE_YX) + + assert triangle.contains(np.array([point, swapped])).tolist() == [True, False] + + +def test_triangle_contains_mask_is_not_reflected(): + """ + Regression: `triangle_contains_mask` unpacked its own vertices as `y, x` + (element 1 as "x") while testing a triangle centroid built as `x, y` + (element 0 as "x"), so it tested the reflected triangle. It reported the + genuinely-interior point (2.0, 0.2) as outside and its reflection + (0.2, 2.0) as inside. + """ + triangle = Triangle(*ASYMMETRIC_TRIANGLE_YX) + + assert triangle.triangle_contains_mask( + triangle_array_with_centroid(POINT_INSIDE_YX) + ) == np.array([True]) + assert triangle.triangle_contains_mask( + triangle_array_with_centroid(POINT_INSIDE_YX[::-1]) + ) == np.array([False]) + + +CIRCLE = Circle(x=1.0, y=2.0, radius=0.5) +TRIANGLE = Triangle(*ASYMMETRIC_TRIANGLE_YX) +POLYGON = Polygon([(0.0, 0.0), (2.0, 0.0), (2.0, 1.0), (0.0, 1.0)]) +SQUARE = Square(top=0.0, bottom=1.0, left=0.0, right=2.0) + + +@pytest.mark.parametrize( + "shape, inside, outside", + [ + ( + CIRCLE, + [(1.0, 2.0), (1.4, 2.0), (1.0, 2.5), (1.3, 2.3)], + [(1.6, 2.0), (1.0, 2.6), (2.0, 3.0)], + ), + ( + TRIANGLE, + [(0.1, 0.1), (2.0, 0.2), (0.0, 0.0), (3.0, 0.0)], + [(2.0, 0.6), (-0.1, 0.5), (4.0, 0.0), (0.0, 1.5)], + ), + ( + POLYGON, + [(1.0, 0.5), (0.0, 0.0), (2.0, 1.0)], + [(3.0, 0.5), (1.0, 1.5), (-0.5, 0.5)], + ), + ( + SQUARE, + [(1.0, 0.5), (0.0, 0.0), (2.0, 1.0)], + [(2.5, 0.5), (1.0, 1.5), (-0.5, 0.5)], + ), + ], +) +def test_contains(shape, inside, outside): + assert shape.contains(np.array(inside)).all() + assert not shape.contains(np.array(outside)).any() + + +def test_square_contains_is_robust_to_coordinate_ordering(): + """ + `mask` and `area` assume `top < bottom` numerically (coordinates from the + top-left corner of the image). `contains` and `boundary` sort the pair, so + they behave identically for a square built from arcsec `(y, x)` + coordinates, where the top of the image is the larger first coordinate. + """ + flipped = Square(top=1.0, bottom=0.0, left=2.0, right=0.0) + + points = np.array([(1.0, 0.5), (0.0, 0.0), (2.5, 0.5), (1.0, 1.5)]) + + assert flipped.contains(points).tolist() == SQUARE.contains(points).tolist() + assert flipped.boundary().tolist() == SQUARE.boundary().tolist() + + +@pytest.mark.parametrize("shape", [CIRCLE, TRIANGLE, POLYGON, SQUARE]) +def test_boundary_is_closed_and_its_edges_lie_inside(shape): + boundary = shape.boundary() + + assert boundary.ndim == 2 + assert boundary.shape[1] == 2 + assert boundary[0] == pytest.approx(boundary[-1]) + + midpoints = 0.5 * (boundary[:-1] + boundary[1:]) + + assert shape.contains(midpoints).all() + + +def test_circle_boundary_sample_count(): + assert CIRCLE.boundary(n=8).shape == (9, 2) + assert CIRCLE.boundary().shape == (101, 2) + + +@pytest.mark.parametrize( + "point", + [ + (1.0, 2.0), + (1.4, 2.0), + (1.6, 2.0), + (1.0, 2.6), + (0.0, 0.0), + ], +) +def test_circle_contains_agrees_with_mask(point): + """ + A triangle small enough to sit at a single coordinate is kept by + `Circle.mask` exactly when that coordinate is inside the circle. + """ + assert ( + CIRCLE.contains(np.array([point]))[0] + == CIRCLE.mask(triangle_array_with_centroid(point))[0] + ) + + +def test_point_contains_raises(): + point = Point(1.0, 2.0) + + with pytest.raises(NotImplementedError): + point.contains(np.array([[1.0, 2.0]])) + + +def test_point_boundary_is_a_single_row(): + assert Point(1.0, 2.0).boundary().tolist() == [[1.0, 2.0]] From 67b1cb754136a2b2717fb554b03f7244a6cdaefc Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 2 Sep 2026 18:08:19 -0400 Subject: [PATCH 2/3] fix(triangles): tile the requested rectangle, not its transpose (image-source-mappings phase 2a) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `CoordinateArrayTriangles{,Np}.for_limits_and_scale` named its first pair of limits `x_min`/`x_max` and tiled them along element 0 of every vertex, while element 0 is the `y` coordinate everywhere else — `ArrayTrianglesNp.for_limits_and_scale`, every PyAuto `(y, x)` grid, and the `element 0 <-> element 0` convention `Shape.contains` / `Shape.mask` are documented with. Both the keyword caller (`AbstractSolver._initial_triangles`) and the positional caller (`AbstractTriangles.for_grid`) therefore tiled the transposed rectangle. A square grid hides this, which is why it shipped; on a 24x80 grid of 0.05" pixels `PointSolver.solve` found one of an Isothermal's two images instead of both, the missing one lying well inside the grid. Also removes a stray `bbbb` statement (a `NameError` on every call) left in `ArrayTrianglesNp.with_vertices` since 2025-11-13, uncaught because the solver only ever reaches the `CoordinateArrayTrianglesNp` override of that name. Found by the ShapeSolver validation suite of image-source-mappings phase 2. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_011EVnQ4sEYrRM7GCvFyveUA --- autoarray/structures/triangles/array_np.py | 1 - .../structures/triangles/coordinate_array.py | 37 +++++++++-- .../triangles/coordinate_array_np.py | 37 +++++++++-- .../structures/triangles/test_coordinate.py | 64 +++++++++++++++++++ 4 files changed, 124 insertions(+), 15 deletions(-) diff --git a/autoarray/structures/triangles/array_np.py b/autoarray/structures/triangles/array_np.py index 4606ee69f..9a623294d 100644 --- a/autoarray/structures/triangles/array_np.py +++ b/autoarray/structures/triangles/array_np.py @@ -207,7 +207,6 @@ def with_vertices(self, vertices: np.ndarray) -> "ArrayTrianglesNp": ------- The new set of triangles with the new vertices. """ - bbbb return ArrayTrianglesNp( indices=self.indices, vertices=vertices, diff --git a/autoarray/structures/triangles/coordinate_array.py b/autoarray/structures/triangles/coordinate_array.py index 2d674c58e..1fd70fbd0 100644 --- a/autoarray/structures/triangles/coordinate_array.py +++ b/autoarray/structures/triangles/coordinate_array.py @@ -46,23 +46,46 @@ def __init__( @classmethod def for_limits_and_scale( cls, - x_min: float, - x_max: float, y_min: float, y_max: float, + x_min: float, + x_max: float, scale: float = 1.0, **_, ): + """ + Tile the rectangle ``[y_min, y_max] x [x_min, x_max]`` with equilateral triangles. + + Element ``0`` of every vertex spans the ``y`` limits and element ``1`` the ``x`` + limits, matching `ArrayTrianglesNp.for_limits_and_scale`, the ``(y, x)`` order of + every PyAuto grid, and the ``element 0 <-> element 0`` convention `Shape.contains` + and `Shape.mask` are documented with. + + The two axes were previously the other way round while the signature named its + first pair ``x_min``/``x_max``, so both the keyword callers (`AbstractSolver`) and + the positional caller (`AbstractTriangles.for_grid`) tiled the *transposed* + rectangle. On a square grid that is invisible; on a rectangular one the solver + searched a box the data does not occupy and silently missed multiple images which + lay inside the grid (a 24x80 grid of 0.05" pixels found one of an Isothermal's two + images instead of both). + + Parameters + ---------- + y_min, y_max, x_min, x_max + The limits of the rectangle to tile. + scale + The side length of the triangles. + """ import jax.numpy as jnp - x_shift = int(2 * x_min / scale) - y_shift = int(y_min / (HEIGHT_FACTOR * scale)) + y_shift = int(2 * y_min / scale) + x_shift = int(x_min / (HEIGHT_FACTOR * scale)) coordinates = [] - for x in range(x_shift, int(2 * x_max / scale) + 1): - for y in range(y_shift - 1, int(y_max / (HEIGHT_FACTOR * scale)) + 2): - coordinates.append([x, y]) + for y in range(y_shift, int(2 * y_max / scale) + 1): + for x in range(x_shift - 1, int(x_max / (HEIGHT_FACTOR * scale)) + 2): + coordinates.append([y, x]) return cls( coordinates=jnp.array(coordinates), diff --git a/autoarray/structures/triangles/coordinate_array_np.py b/autoarray/structures/triangles/coordinate_array_np.py index 21b9a0dc1..511cd7e9e 100644 --- a/autoarray/structures/triangles/coordinate_array_np.py +++ b/autoarray/structures/triangles/coordinate_array_np.py @@ -92,21 +92,44 @@ def flip_array(self) -> np.ndarray: @classmethod def for_limits_and_scale( cls, - x_min: float, - x_max: float, y_min: float, y_max: float, + x_min: float, + x_max: float, scale: float = 1.0, **_, ): - x_shift = int(2 * x_min / scale) - y_shift = int(y_min / (HEIGHT_FACTOR * scale)) + """ + Tile the rectangle ``[y_min, y_max] x [x_min, x_max]`` with equilateral triangles. + + Element ``0`` of every vertex spans the ``y`` limits and element ``1`` the ``x`` + limits, matching `ArrayTrianglesNp.for_limits_and_scale`, the ``(y, x)`` order of + every PyAuto grid, and the ``element 0 <-> element 0`` convention `Shape.contains` + and `Shape.mask` are documented with. + + The two axes were previously the other way round while the signature named its + first pair ``x_min``/``x_max``, so both the keyword callers (`AbstractSolver`) and + the positional caller (`AbstractTriangles.for_grid`) tiled the *transposed* + rectangle. On a square grid that is invisible; on a rectangular one the solver + searched a box the data does not occupy and silently missed multiple images which + lay inside the grid (a 24x80 grid of 0.05" pixels found one of an Isothermal's two + images instead of both). + + Parameters + ---------- + y_min, y_max, x_min, x_max + The limits of the rectangle to tile. + scale + The side length of the triangles. + """ + y_shift = int(2 * y_min / scale) + x_shift = int(x_min / (HEIGHT_FACTOR * scale)) coordinates = [] - for x in range(x_shift, int(2 * x_max / scale) + 1): - for y in range(y_shift - 1, int(y_max / (HEIGHT_FACTOR * scale)) + 2): - coordinates.append([x, y]) + for y in range(y_shift, int(2 * y_max / scale) + 1): + for x in range(x_shift - 1, int(x_max / (HEIGHT_FACTOR * scale)) + 2): + coordinates.append([y, x]) return CoordinateArrayTrianglesNp( coordinates=np.array(coordinates, dtype=np.int32), diff --git a/test_autoarray/structures/triangles/test_coordinate.py b/test_autoarray/structures/triangles/test_coordinate.py index 84a70434f..1a60ed640 100644 --- a/test_autoarray/structures/triangles/test_coordinate.py +++ b/test_autoarray/structures/triangles/test_coordinate.py @@ -325,3 +325,67 @@ def test_area(one_triangle): assert neighborhood.area == 4 * ONE_TRIANGLE_AREA assert neighborhood.up_sample().area == 4 * ONE_TRIANGLE_AREA assert neighborhood.neighborhood().area == 10 * ONE_TRIANGLE_AREA + + +def test_for_limits_and_scale__element_0_spans_the_y_limits(): + """ + Regression: `for_limits_and_scale` named its first pair of limits ``x_min``/``x_max`` + and tiled them along element ``0`` of every vertex, while element ``0`` is the ``y`` + coordinate everywhere else (``ArrayTrianglesNp``, every PyAuto ``(y, x)`` grid, and + the ``element 0 <-> element 0`` convention of ``Shape.contains`` / ``Shape.mask``). + A rectangular extent was therefore tiled transposed, so the point solver searched a + box the data does not occupy. A square extent hides it, which is why it shipped. + """ + triangles = CoordinateArrayTrianglesNp.for_limits_and_scale( + y_min=-2.0, + y_max=2.0, + x_min=-1.0, + x_max=1.0, + scale=0.5, + ) + + vertices = triangles.vertices + + assert min(vertices[:, 0]) <= -2.0 + assert max(vertices[:, 0]) >= 2.0 + + assert min(vertices[:, 1]) <= -1.0 + assert max(vertices[:, 1]) >= 1.0 + + # The transposed tiling would have run element 1 out to +/- 2.0. + assert max(vertices[:, 1]) < 2.0 + + +def test_for_grid__rectangular_grid_is_tiled_the_same_way_round(): + """ + ``AbstractTriangles.for_grid`` passes the limits *positionally* as + ``(y_min, y_max, x_min, x_max)``, so the signature order matters as much as the + tiling does. + """ + from autoarray.structures.grids.uniform_2d import Grid2D + + grid = Grid2D.uniform(shape_native=(24, 80), pixel_scales=0.05) + + triangles = CoordinateArrayTrianglesNp.for_grid(grid=grid) + + vertices = triangles.vertices + + assert max(vertices[:, 0]) < max(vertices[:, 1]) + + +def test_array_triangles_with_vertices_returns_new_triangles(): + """ + Regression: ``ArrayTrianglesNp.with_vertices`` carried a stray ``bbbb`` statement (a + ``NameError`` on every call) from 2025-11, uncaught because the solver only ever + reached the ``CoordinateArrayTrianglesNp`` override of the same name. + """ + from autoarray.structures.triangles.array_np import ArrayTrianglesNp + + triangles = ArrayTrianglesNp( + indices=np.array([[0, 1, 2]]), + vertices=np.array([[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]), + ) + + new_vertices = np.array([[0.0, 0.0], [2.0, 0.0], [0.0, 2.0]]) + + assert (triangles.with_vertices(new_vertices).vertices == new_vertices).all() From a9a9120a09719d1990c6f57651d8b21a3c15ab87 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 2 Sep 2026 18:38:24 -0400 Subject: [PATCH 3/3] fix(plot): label every polygon of a multi-image region (image-source-mappings phase 2a) `plot_regions` drew one label per region entry, at the mean of all of that region's polygons. A region whose polygons are the multiple images of one lensed source has them on opposite sides of the lens, so the label landed at their midpoint -- on empty sky, labelling nothing, which is exactly the case the mappings overlay exists for. The label is now repeated once per polygon, at that polygon's own centre. A single-polygon region is unchanged: its polygon's mean is the region's mean. Found by the phase 2 verification renders. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_011EVnQ4sEYrRM7GCvFyveUA --- autoarray/plot/utils.py | 38 +++++++++++++-------- test_autoarray/plot/test_utils.py | 55 +++++++++++++++++++++++++++++++ 2 files changed, 79 insertions(+), 14 deletions(-) diff --git a/autoarray/plot/utils.py b/autoarray/plot/utils.py index 3f8186d04..9587517e3 100644 --- a/autoarray/plot/utils.py +++ b/autoarray/plot/utils.py @@ -1205,8 +1205,10 @@ def plot_regions( region_alpha The alpha of each region's fill; the outline is always drawn opaque. region_labels - A text label drawn at the centre of each region, e.g. ``["1", "2"]``. ``None`` draws no - labels. + A text label drawn at the centre of each polygon of each region, e.g. ``["1", "2"]``. + A region with several polygons -- the multiple images of one lensed source -- gets its + label repeated once per image, so every image is labelled and no label lands on the + empty sky between them. ``None`` draws no labels. """ if regions is None: return @@ -1240,15 +1242,23 @@ def plot_regions( ax.plot(polygon[:, 1], polygon[:, 0], color=color, linewidth=1, zorder=5) if region_labels is not None and i < len(region_labels) and len(points) > 0: - stacked = np.concatenate(points, axis=0) - - ax.annotate( - str(region_labels[i]), - xy=(float(np.mean(stacked[:, 1])), float(np.mean(stacked[:, 0]))), - color=color, - fontsize=12, - fontweight="bold", - ha="center", - va="center", - zorder=6, - ) + # One label per polygon, at that polygon's own centre -- not one label at the + # centre of all of them. A region whose polygons are the multiple images of a + # lensed source has its polygons on opposite sides of the lens, so a single + # label at their combined mean lands between them, on empty sky, labelling + # nothing. A single-polygon region is unaffected: its polygon's mean is the + # region's mean. + for polygon in points: + ax.annotate( + str(region_labels[i]), + xy=( + float(np.mean(polygon[:, 1])), + float(np.mean(polygon[:, 0])), + ), + color=color, + fontsize=12, + fontweight="bold", + ha="center", + va="center", + zorder=6, + ) diff --git a/test_autoarray/plot/test_utils.py b/test_autoarray/plot/test_utils.py index 5d4334cda..2196eed8e 100644 --- a/test_autoarray/plot/test_utils.py +++ b/test_autoarray/plot/test_utils.py @@ -505,3 +505,58 @@ def test__plot_inversion_reconstruction__regions_overlay_is_output( ) assert str(Path(plot_path) / "reconstruction_regions.png") in plot_patch.paths + + +def test__plot_regions__multi_polygon_region_is_labelled_once_per_polygon(): + """ + A region whose polygons are the multiple images of one lensed source has its polygons + on opposite sides of the lens. A single label at their combined mean therefore landed + between the images, on empty sky, labelling nothing -- so the label is repeated once + per polygon, at that polygon's own centre. + """ + import matplotlib.pyplot as plt + + left = np.array([[1.0, -2.0], [1.0, -1.0], [0.0, -1.0], [1.0, -2.0]]) + right = np.array([[1.0, 1.0], [1.0, 2.0], [0.0, 2.0], [1.0, 1.0]]) + + figure, ax = plt.subplots() + + try: + plot_utils.plot_regions(ax=ax, regions=[[left, right]], region_labels=["1"]) + + texts = [text for text in ax.texts if text.get_text() == "1"] + + assert len(texts) == 2 + + x_positions = sorted(float(text.get_position()[0]) for text in texts) + + # One label over each image, not one between them at x = 0. + assert x_positions[0] == pytest.approx(np.mean(left[:, 1])) + assert x_positions[1] == pytest.approx(np.mean(right[:, 1])) + finally: + plt.close(figure) + + +def test__plot_regions__single_polygon_region_is_labelled_once(): + """ + The single-polygon case is unchanged: one label, at the polygon's centre. + """ + import matplotlib.pyplot as plt + + polygon = np.array( + [[1.0, -1.0], [1.0, 1.0], [-1.0, 1.0], [-1.0, -1.0], [1.0, -1.0]] + ) + + figure, ax = plt.subplots() + + try: + plot_utils.plot_regions(ax=ax, regions=[[polygon]], region_labels=["1"]) + + texts = [text for text in ax.texts if text.get_text() == "1"] + + assert len(texts) == 1 + assert float(texts[0].get_position()[0]) == pytest.approx( + np.mean(polygon[:, 1]) + ) + finally: + plt.close(figure)