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/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/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/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) 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() 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]]