From 36b75e0722b8d07441069385ad692a61f9e03309 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Wed, 2 Sep 2026 13:50:43 -0400 Subject: [PATCH] =?UTF-8?q?feat:=20image=E2=86=94source=20mapping=20region?= =?UTF-8?q?s,=20source=20clumps,=20restored=20subplot=5Fmappings?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add autoarray/inversion/mappings/ (Mapping, ImageRegion, mesh-graph connected components, image-plane regions from the mapping matrix, arcsec boundary polygons), Inversion.source_clumps_from / mappings_from, Mapper.mappings_from, a regions= polygon overlay on plot_array / plot_inversion_reconstruction / plot_mapper, and rewrite subplot_mappings as the 2x2 colour-matched figure of how each bright source clump maps to its multiple images. Config keys total_mappings / mappings_threshold / mappings_min_pixels replace the dead total_mappings_pixels (read as a fallback). Phase 1 of the image-source-mappings epic. Closes #515. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01V2ujRnmqc3jGh82hZAQ9iJ --- autoarray/__init__.py | 2 + autoarray/config/visualize/general.yaml | 4 +- autoarray/inversion/inversion/abstract.py | 118 +++ autoarray/inversion/mappers/abstract.py | 46 +- autoarray/inversion/mappings/__init__.py | 7 + autoarray/inversion/mappings/mapping.py | 754 ++++++++++++++++++ autoarray/inversion/plot/inversion_plots.py | 132 ++- autoarray/inversion/plot/mapper_plots.py | 50 +- autoarray/plot/array.py | 42 +- autoarray/plot/inversion.py | 63 +- autoarray/plot/utils.py | 120 ++- .../inversion/inversion/test_abstract.py | 107 +++ test_autoarray/inversion/mappings/__init__.py | 0 .../inversion/mappings/test_mapping.py | 237 ++++++ .../inversion/plot/test_inversion_plotters.py | 34 + .../inversion/plot/test_mapper_plotters.py | 21 + test_autoarray/plot/test_utils.py | 36 +- 17 files changed, 1732 insertions(+), 41 deletions(-) create mode 100644 autoarray/inversion/mappings/__init__.py create mode 100644 autoarray/inversion/mappings/mapping.py create mode 100644 test_autoarray/inversion/mappings/__init__.py create mode 100644 test_autoarray/inversion/mappings/test_mapping.py diff --git a/autoarray/__init__.py b/autoarray/__init__.py index 314e5cc7a..c96ab9ea1 100644 --- a/autoarray/__init__.py +++ b/autoarray/__init__.py @@ -37,6 +37,8 @@ from .inversion.mesh.border_relocator import BorderRelocator from .inversion.pixelization import Pixelization from .inversion.mappers.abstract import Mapper +from .inversion.mappings.mapping import ImageRegion +from .inversion.mappings.mapping import Mapping from .inversion.mesh.image_mesh.abstract import AbstractImageMesh from .inversion.mesh.mesh.abstract import AbstractMesh from .inversion.mesh.interpolator.rectangular import InterpolatorRectangular diff --git a/autoarray/config/visualize/general.yaml b/autoarray/config/visualize/general.yaml index 7098109c1..04e1d1fcb 100644 --- a/autoarray/config/visualize/general.yaml +++ b/autoarray/config/visualize/general.yaml @@ -8,7 +8,9 @@ general: output_format: show # Default output format: "show" displays the figure interactively via plt.show(), "png"/"pdf"/etc. saves to file. inversion: reconstruction_vmax_factor: 0.5 - total_mappings_pixels: 8 # The number of source pixels used when plotting the subplot_mappings of a pixelization. + total_mappings: 5 # The maximum number of source clumps drawn by subplot_mappings. + mappings_threshold: 0.5 # A source pixel joins a clump if its reconstructed value exceeds this fraction of the reconstruction's maximum. + mappings_min_pixels: 3 # Connected groups of source pixels smaller than this are not drawn as a clump. zoom: plane_percent: 0.01 inversion_percent: 0.01 # Plots of an Inversion's reconstruction use the reconstructed data's bright value multiplied by this factor. diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py index 43b915100..68b9be306 100644 --- a/autoarray/inversion/inversion/abstract.py +++ b/autoarray/inversion/inversion/abstract.py @@ -1311,3 +1311,121 @@ def max_pixel_centre(self, mapper_index: int = 0) -> Grid2DIrregular: ) return max_pixel_centre + + def source_clumps_from( + self, + mapper_index: int = 0, + threshold: float = 0.5, + min_pixels: int = 3, + total_clumps: Optional[int] = None, + pix_indexes: Optional[List] = None, + ) -> List[np.ndarray]: + """ + Returns the clumps of the reconstruction, each a group of connected mesh pixels which are + bright relative to the reconstruction's maximum value. + + Mesh pixels whose reconstructed value exceeds `threshold * max(reconstruction)` are kept, + split into connected components over the mesh neighbour graph, filtered by `min_pixels` and + ordered brightest first. + + In gravitational lensing a clump is a distinct bright structure of the source galaxy. The + `threshold` is therefore the knob which decides what counts as a distinct structure: ~0.5 + isolates one smooth source, ~0.2 merges two nearby galaxies into one clump and ~0.8 splits + a source into its individual star-forming knots. + + Parameters + ---------- + mapper_index + The index of the mapper in the inversion whose reconstruction is clumped, where there + may be multiple mappers in the inversion. + threshold + Mesh pixels are in a clump if their reconstructed value exceeds this fraction of the + reconstruction's maximum value. + min_pixels + Connected groups with fewer than this many mesh pixels are discarded. + total_clumps + The maximum number of clumps returned, keeping the brightest. `None` returns all of them. + pix_indexes + An explicit list of mesh pixel index groups which bypasses the clump finding entirely, + returning them unchanged. + + Returns + ------- + A list of 1D integer arrays, one per clump, ordered by decreasing peak reconstructed value. + """ + from autoarray.inversion.mappings.mapping import connected_components_from + from autoarray.inversion.mappings.mapping import pix_index_groups_from + + if pix_indexes is not None: + return pix_index_groups_from(pix_indexes=pix_indexes) + + mapper = self.cls_list_from(cls=Mapper)[mapper_index] + + reconstruction = np.asarray(self.reconstruction_dict[mapper]) + + if reconstruction.size == 0: + return [] + + max_value = float(np.max(reconstruction)) + + if max_value <= 0.0: + return [] + + indexes = np.where(reconstruction > threshold * max_value)[0] + + if indexes.size == 0: + return [] + + clumps = connected_components_from(indexes=indexes, neighbors=mapper.neighbors) + + clumps = [clump for clump in clumps if clump.shape[0] >= min_pixels] + + clumps.sort(key=lambda clump: -float(np.max(reconstruction[clump]))) + + if total_clumps is not None: + clumps = clumps[:total_clumps] + + return clumps + + def mappings_from( + self, + mapper_index: int = 0, + weight_threshold: float = 0.0, + **clump_kwargs, + ) -> List["Mapping"]: + """ + Returns the `Mapping` of every clump of the reconstruction, pairing each source-plane clump + with the image-plane regions (the multiple images) it maps to. + + Parameters + ---------- + mapper_index + The index of the mapper in the inversion whose reconstruction is mapped, where there may + be multiple mappers in the inversion. + weight_threshold + A data pixel is in an image-plane region if its summed mapping weight to the clump + exceeds this value. + clump_kwargs + The `threshold`, `min_pixels`, `total_clumps` and `pix_indexes` inputs + of `source_clumps_from`. + + Returns + ------- + One `Mapping` per source-plane clump, ordered by decreasing peak reconstructed value. + """ + from autoarray.inversion.mappings import mapping as mapping_module + + mapper = self.cls_list_from(cls=Mapper)[mapper_index] + + reconstruction = np.asarray(self.reconstruction_dict[mapper]) + + clumps = self.source_clumps_from(mapper_index=mapper_index, **clump_kwargs) + + peak_values = [float(np.max(reconstruction[clump])) for clump in clumps] + + return mapping_module.mappings_from( + mapper=mapper, + pix_indexes=clumps, + weight_threshold=weight_threshold, + peak_values=peak_values, + ) diff --git a/autoarray/inversion/mappers/abstract.py b/autoarray/inversion/mappers/abstract.py index 2259213b3..45bbcd602 100644 --- a/autoarray/inversion/mappers/abstract.py +++ b/autoarray/inversion/mappers/abstract.py @@ -619,7 +619,9 @@ def extent_from( x_min, x_max, y_min, y_max = extent x_centre = 0.5 * (x_min + x_max) y_centre = 0.5 * (y_min + y_max) - target_half = 0.5 * max(x_max - x_min, y_max - y_min) * zoom_extent_scale + target_half = ( + 0.5 * max(x_max - x_min, y_max - y_min) * zoom_extent_scale + ) bound = geometry_util.extent_symmetric_from( extent=self.source_plane_mesh_grid.geometry.extent @@ -630,8 +632,8 @@ def extent_from( y_centre - bound[2], bound[3] - y_centre, ) - bound_cap_half = 0.7 * 0.5 * min( - bound[1] - bound[0], bound[3] - bound[2] + bound_cap_half = ( + 0.7 * 0.5 * min(bound[1] - bound[0], bound[3] - bound[2]) ) final_half = min(target_half, max_allowable_half, bound_cap_half) @@ -668,3 +670,41 @@ def mesh_pixels_per_image_pixels(self): values=mesh_pixels_per_image_pixels, mask=self.mask, ) + + def mappings_from( + self, + pix_indexes, + weight_threshold: float = 0.0, + min_pixels: int = 1, + ) -> List["Mapping"]: + """ + Returns the `Mapping` of one or more groups of mesh pixels, pairing each source-plane group + with the image-plane regions (the multiple images) it maps to. + + This is the bare-`Mapper` route into the mapping objects, used when there is no inversion to + find clumps in and the mesh pixels of interest are chosen by hand (e.g. in a tutorial). Each + returned `Mapping` therefore has a `peak_value` of `None`. + + Parameters + ---------- + pix_indexes + The mesh pixel indexes, either a flat sequence (one group) or a nested sequence (one + group per entry). + weight_threshold + A data pixel is in an image-plane region if its summed mapping weight to the group + exceeds this value. + min_pixels + Connected image-plane regions with fewer than this many pixels are discarded. + + Returns + ------- + One `Mapping` per group of mesh pixels. + """ + from autoarray.inversion.mappings import mapping as mapping_module + + return mapping_module.mappings_from( + mapper=self, + pix_indexes=pix_indexes, + weight_threshold=weight_threshold, + min_pixels=min_pixels, + ) diff --git a/autoarray/inversion/mappings/__init__.py b/autoarray/inversion/mappings/__init__.py new file mode 100644 index 000000000..e7f58cbdf --- /dev/null +++ b/autoarray/inversion/mappings/__init__.py @@ -0,0 +1,7 @@ +from autoarray.inversion.mappings.mapping import ImageRegion +from autoarray.inversion.mappings.mapping import Mapping +from autoarray.inversion.mappings.mapping import connected_components_from +from autoarray.inversion.mappings.mapping import contours_from_bool_native +from autoarray.inversion.mappings.mapping import image_regions_from +from autoarray.inversion.mappings.mapping import image_regions_from_slim_mask +from autoarray.inversion.mappings.mapping import source_contours_from diff --git a/autoarray/inversion/mappings/mapping.py b/autoarray/inversion/mappings/mapping.py new file mode 100644 index 000000000..478266028 --- /dev/null +++ b/autoarray/inversion/mappings/mapping.py @@ -0,0 +1,754 @@ +""" +Plotting-agnostic objects describing how a region of the source-plane maps to the image-plane. + +A `Mapping` pairs one source-plane region (a group of mesh pixels, e.g. a clump of the +reconstruction) with the list of image-plane connected regions (the multiple images) that the +region maps to. Both planes are described as polygons in scaled (e.g. arc-second) coordinates, +so that they can be overlaid on any plot without depending on how the underlying data is binned, +zoomed or cropped. + +Everything in this module is pure numpy and is never jitted, because it is diagnostic / +visualization code which runs once per figure. +""" + +from dataclasses import dataclass +from typing import List, Optional, Sequence, Tuple + +import numpy as np + + +COLORS_DEFAULT = ["r", "g", "b", "m", "c", "y"] + + +@dataclass(frozen=True, eq=False) +class ImageRegion: + """ + A single connected region of image-plane pixels which a source-plane region maps to. + + For a strong lens this is one of the multiple images of a source clump. + + Parameters + ---------- + slim_indexes + The slim indexes (indexes into the masked data, without over-sampling) of every image + pixel in the region. + mask + A `Mask2D` whose `False` entries are the pixels inside the region, using the same 2D + shape, pixel scales and origin as the data the region was computed from. + contours + The boundary loops of the region, each an ``(N, 2)`` array of ``(y, x)`` scaled + coordinates whose first and last entries are identical (a closed loop). A region with a + hole, or which pinches at a corner, has more than one loop. + centre + The ``(y, x)`` scaled coordinate of the geometric centre (unweighted mean) of the region. + """ + + slim_indexes: np.ndarray + mask: "Mask2D" + contours: List[np.ndarray] + centre: Tuple[float, float] + + @property + def pixel_coordinates(self) -> np.ndarray: + """ + The ``(y, x)`` native pixel coordinates of every pixel inside the region, shape ``(N, 2)``. + """ + return np.argwhere(~np.asarray(self.mask)) + + @property + def scaled_coordinates(self) -> np.ndarray: + """ + The ``(y, x)`` scaled coordinates of the centre of every pixel inside the region, shape + ``(N, 2)``. + """ + geometry = self.mask.geometry + + pixels = self.pixel_coordinates + + y = geometry.scaled_maxima[0] - (pixels[:, 0] + 0.5) * geometry.pixel_scales[0] + x = geometry.scaled_minima[1] + (pixels[:, 1] + 0.5) * geometry.pixel_scales[1] + + return np.stack([y, x], axis=-1) + + def area(self) -> float: + """ + The area of the region in scaled units squared (e.g. arcsec^2), computed as its number of + pixels multiplied by the area of a single pixel. + + Returns + ------- + The scaled area of the region. + """ + pixel_scales = self.mask.geometry.pixel_scales + + return float(len(self.slim_indexes) * pixel_scales[0] * pixel_scales[1]) + + def values_from(self, array) -> np.ndarray: + """ + The values of `array` at every pixel inside the region. + + Parameters + ---------- + array + An `Array2D` (or 2D numpy array in native format) defined on the same mask geometry as + the region. + + Returns + ------- + The values of `array` inside the region, ordered as `pixel_coordinates`. + """ + native = ( + np.asarray(array.native.array) + if hasattr(array, "native") + else np.asarray(array) + ) + + if native.ndim != 2: + raise ValueError( + "The array input to an ImageRegion must be an Array2D or a 2D native numpy array, " + f"but an array of shape {native.shape} was input." + ) + + pixels = self.pixel_coordinates + + return native[pixels[:, 0], pixels[:, 1]] + + def brightest_coordinate_from(self, array) -> Tuple[float, float]: + """ + The ``(y, x)`` scaled coordinate of the brightest pixel of `array` inside the region. + + Parameters + ---------- + array + An `Array2D` (or 2D native numpy array) defined on the same mask geometry as the region. + + Returns + ------- + The ``(y, x)`` scaled coordinate of the region's brightest pixel. + """ + values = self.values_from(array=array) + + index = int(np.argmax(values)) + + coordinate = self.scaled_coordinates[index] + + return float(coordinate[0]), float(coordinate[1]) + + def centroid_from(self, array) -> Tuple[float, float]: + """ + The flux-weighted ``(y, x)`` scaled centroid of `array` inside the region. + + Parameters + ---------- + array + An `Array2D` (or 2D native numpy array) defined on the same mask geometry as the region. + + Returns + ------- + The ``(y, x)`` scaled coordinate of the region's flux-weighted centroid. + """ + values = self.values_from(array=array) + + total = float(np.sum(values)) + + if total == 0.0: + return self.centre + + coordinates = self.scaled_coordinates + + y = float(np.sum(coordinates[:, 0] * values) / total) + x = float(np.sum(coordinates[:, 1] * values) / total) + + return y, x + + def flux_from(self, array) -> float: + """ + The summed value of `array` over every pixel inside the region. + + Parameters + ---------- + array + An `Array2D` (or 2D native numpy array) defined on the same mask geometry as the region. + + Returns + ------- + The total flux of the region. + """ + return float(np.sum(self.values_from(array=array))) + + +@dataclass(frozen=True, eq=False) +class Mapping: + """ + One source-plane region and the image-plane regions (multiple images) it maps to. + + Parameters + ---------- + pix_indexes + The indexes of the mesh pixels in the source-plane region. + source_contours + The boundaries of the source-plane region, each an ``(N, 2)`` array of ``(y, x)`` scaled + coordinates forming a closed loop. One loop per mesh pixel (rectangular cells / Voronoi + cells), so the region is drawn as the union of its pixels. + source_centre + The ``(y, x)`` scaled coordinate of the mean of the region's mesh pixel centres. + image_regions + The connected image-plane regions the source region maps to. + peak_value + The maximum reconstructed value over the region's mesh pixels, or `None` when the mapping + was computed from a bare `Mapper` (which has no reconstruction). + """ + + pix_indexes: np.ndarray + source_contours: List[np.ndarray] + source_centre: Tuple[float, float] + image_regions: List[ImageRegion] + peak_value: Optional[float] = None + + @property + def image_contours(self) -> List[np.ndarray]: + """ + Every image-plane boundary loop of the mapping, flattened across its image regions. + """ + return [contour for region in self.image_regions for contour in region.contours] + + +def pix_index_groups_from(pix_indexes: Sequence) -> List[np.ndarray]: + """ + Normalise a user-input `pix_indexes` into a list of mesh-pixel index groups. + + A flat sequence of indexes (e.g. ``[0, 1, 2]``) is a single group; a nested sequence + (e.g. ``[[0, 1], [5]]``) is one group per entry. + + Parameters + ---------- + pix_indexes + The mesh pixel indexes, either flat or nested. + + Returns + ------- + A list of 1D integer arrays, one per group. + """ + groups = list(pix_indexes) + + if len(groups) == 0: + return [] + + if np.ndim(groups[0]) == 0: + return [np.asarray(groups, dtype=int)] + + return [np.asarray(group, dtype=int) for group in groups] + + +def connected_components_from(indexes, neighbors) -> List[np.ndarray]: + """ + Split `indexes` into connected components of the mesh neighbour graph. + + Two mesh pixels are in the same component if a path of neighbouring pixels between them exists + which passes only through pixels in `indexes`. + + Parameters + ---------- + indexes + The mesh pixel indexes to group. + neighbors + A `Neighbors` object mapping every mesh pixel to the indexes of its neighbours, whose + `sizes` attribute gives the number of valid entries of each row. + + Returns + ------- + A list of 1D integer arrays, one per connected component, each sorted ascending. Components + are ordered by the first appearance of one of their members in `indexes`. + """ + indexes = np.asarray(indexes, dtype=int).ravel() + + neighbors_arr = np.asarray(neighbors) + sizes = np.asarray(neighbors.sizes, dtype=int) + + in_group = set(int(index) for index in indexes) + visited = set() + + components = [] + + for index in indexes: + index = int(index) + + if index in visited: + continue + + visited.add(index) + + stack = [index] + component = [] + + while stack: + pixel = stack.pop() + component.append(pixel) + + row = neighbors_arr[pixel][: sizes[pixel]] + + for neighbor in row: + neighbor = int(neighbor) + + if neighbor < 0 or neighbor in visited or neighbor not in in_group: + continue + + visited.add(neighbor) + stack.append(neighbor) + + components.append(np.sort(np.asarray(component, dtype=int))) + + return components + + +def contours_from_bool_native(bool_native: np.ndarray, geometry) -> List[np.ndarray]: + """ + The boundary loops of the `True` pixels of a 2D boolean array, in scaled coordinates. + + Every edge of a `True` pixel which is not shared with another `True` pixel is a boundary edge. + The boundary edges are directed so that the region is always on the same side of them, and are + then chained end-to-end into closed loops. A region with a hole emits the hole as an extra + loop, and a region which pinches to a point at a corner emits one loop per lobe. + + Parameters + ---------- + bool_native + A 2D boolean array in native format, whose `True` entries are the region. + geometry + A `Geometry2D` (e.g. `mask.geometry`) giving the `pixel_scales`, `scaled_maxima` and + `scaled_minima` used to convert pixel corners to scaled coordinates. + + Returns + ------- + A list of ``(N, 2)`` arrays of ``(y, x)`` scaled coordinates, each a closed loop whose first + and last entries are identical. + """ + bool_native = np.asarray(bool_native, dtype=bool) + + padded = np.zeros( + shape=(bool_native.shape[0] + 2, bool_native.shape[1] + 2), dtype=bool + ) + padded[1:-1, 1:-1] = bool_native + + inside = padded[1:-1, 1:-1] + + # Directed boundary edges in pixel-corner space, each stored as start -> end. The directions + # are chosen so that the region is consistently on one side, which makes the chaining below + # trace each loop without ambiguity except at pinch vertices (where either choice is valid). + + edges = {} + + def _add(start, end): + edges.setdefault(start, []).append(end) + + rows, cols = np.nonzero(inside) + + for row, col in zip(rows.tolist(), cols.tolist()): + if not padded[row, col + 1]: + _add((row, col), (row, col + 1)) + if not padded[row + 1, col + 2]: + _add((row, col + 1), (row + 1, col + 1)) + if not padded[row + 2, col + 1]: + _add((row + 1, col + 1), (row + 1, col)) + if not padded[row + 1, col]: + _add((row + 1, col), (row, col)) + + loops = [] + + starts = list(edges.keys()) + + for start in starts: + while edges.get(start): + loop = [start] + point = start + + while True: + ends = edges.get(point) + + if not ends: + raise ValueError( + "The boundary edges of a region could not be chained into a closed loop, " + "which should be impossible for a boolean pixel region." + ) + + point = ends.pop() + loop.append(point) + + if point == start: + break + + loops.append(np.asarray(loop, dtype=int)) + + y_max = geometry.scaled_maxima[0] + x_min = geometry.scaled_minima[1] + pixel_scales = geometry.pixel_scales + + contours = [] + + for loop in loops: + y = y_max - loop[:, 0] * pixel_scales[0] + x = x_min + loop[:, 1] * pixel_scales[1] + + contours.append(np.stack([y, x], axis=-1)) + + return contours + + +def image_regions_from_slim_mask( + mask, slim_bool, min_pixels: int = 1 +) -> List[ImageRegion]: + """ + The connected image-plane regions of a boolean flag over the masked data pixels. + + The slim boolean is lifted to native 2D, split into connected components with 8-connectivity + (`scipy.ndimage.label`), and each component larger than `min_pixels` is returned as an + `ImageRegion` with its boundary contours in scaled coordinates. + + Parameters + ---------- + mask + The `Mask2D` of the data, whose `False` entries are the (slim ordered) data pixels. + slim_bool + A 1D boolean array of length `mask.pixels_in_mask`, `True` for pixels in the region. + min_pixels + Connected components with fewer than this many pixels are discarded. + + Returns + ------- + The connected regions, ordered largest first. + """ + from scipy import ndimage + + from autoarray.mask.mask_2d import Mask2D + + mask_arr = np.asarray(mask) + + slim_bool = np.asarray(slim_bool, dtype=bool) + + region_native = np.zeros(shape=mask_arr.shape, dtype=bool) + region_native[~mask_arr] = slim_bool + + # The slim index of every unmasked pixel, in the native raster order the slim ordering uses. + slim_index_native = np.full(shape=mask_arr.shape, fill_value=-1, dtype=int) + slim_index_native[~mask_arr] = np.arange(int(np.sum(~mask_arr))) + + labels, total_labels = ndimage.label( + region_native, structure=np.ones(shape=(3, 3), dtype=int) + ) + + regions = [] + + for label in range(1, total_labels + 1): + label_native = labels == label + + slim_indexes = np.sort(slim_index_native[label_native]) + + if slim_indexes.shape[0] < min_pixels: + continue + + region_mask = Mask2D( + mask=~label_native, + pixel_scales=mask.pixel_scales, + origin=mask.origin, + ) + + contours = contours_from_bool_native( + bool_native=label_native, geometry=mask.geometry + ) + + pixels = np.argwhere(label_native) + + geometry = mask.geometry + + centre = ( + float( + np.mean( + geometry.scaled_maxima[0] + - (pixels[:, 0] + 0.5) * geometry.pixel_scales[0] + ) + ), + float( + np.mean( + geometry.scaled_minima[1] + + (pixels[:, 1] + 0.5) * geometry.pixel_scales[1] + ) + ), + ) + + regions.append( + ImageRegion( + slim_indexes=slim_indexes, + mask=region_mask, + contours=contours, + centre=centre, + ) + ) + + regions.sort(key=lambda region: -region.slim_indexes.shape[0]) + + return regions + + +def image_regions_from( + mapper, + pix_indexes, + weight_threshold: float = 0.0, + min_pixels: int = 1, +) -> List[ImageRegion]: + """ + The connected image-plane regions which a group of mesh pixels maps to. + + The mapper's `mapping_matrix` (shape ``[data_pixels, mesh_pixels]``, already folded over + over-sampled sub-pixels) is summed over `pix_indexes`, giving the total weight with which every + data pixel maps to the group. Data pixels whose weight exceeds `weight_threshold` form the + image-plane region, which is then split into connected components. + + Parameters + ---------- + mapper + The `Mapper` whose mappings between the data and the mesh are used. + pix_indexes + The indexes of the mesh pixels forming the source-plane region. + weight_threshold + A data pixel is in the image-plane region if its summed mapping weight exceeds this value. + Raising it keeps only the data pixels which map dominantly to the source region. + min_pixels + Connected image-plane regions with fewer than this many pixels are discarded. + + Returns + ------- + The connected image-plane regions, ordered largest first. + """ + pix_indexes = np.asarray(pix_indexes, dtype=int).ravel() + + mapping_matrix = np.asarray(mapper.mapping_matrix) + + weights = mapping_matrix[:, pix_indexes].sum(axis=1) + + return image_regions_from_slim_mask( + mask=mapper.mask, + slim_bool=weights > weight_threshold, + min_pixels=min_pixels, + ) + + +def _rectangular_edges_from(mapper) -> Tuple[np.ndarray, np.ndarray]: + """ + The ``(y, x)`` cell edge coordinates of a rectangular mesh, one more edge than pixels per axis. + + Uniform rectangular meshes derive their edges from the mesh geometry's pixel scales and origin + (matching the `imshow` extent the reconstruction is drawn with), whereas adaptive rectangular + meshes use the mesh geometry's `edges_transformed` (matching the `pcolormesh` grid). + + Parameters + ---------- + mapper + A `Mapper` whose interpolator is rectangular. + + Returns + ------- + The ``(y_edges, x_edges)`` arrays, of length ``shape[0] + 1`` and ``shape[1] + 1``. + """ + from autoarray.inversion.mesh.interpolator.rectangular_uniform import ( + InterpolatorRectangularUniform, + ) + + mesh_geometry = mapper.mesh_geometry + + shape = mesh_geometry.shape + + if isinstance(mapper.interpolator, InterpolatorRectangularUniform): + geometry = mesh_geometry.geometry + + y_edges = ( + geometry.scaled_maxima[0] + - np.arange(shape[0] + 1) * geometry.pixel_scales[0] + ) + x_edges = ( + geometry.scaled_minima[1] + + np.arange(shape[1] + 1) * geometry.pixel_scales[1] + ) + + return y_edges, x_edges + + y_edges, x_edges = np.asarray(mesh_geometry.edges_transformed).T + + return y_edges, x_edges + + +def _voronoi_contours_from(mapper, pix_indexes: np.ndarray) -> List[np.ndarray]: + """ + The Voronoi cell polygons of a group of Delaunay (or KNN) mesh vertices. + + A Delaunay mesh pixel is a vertex of the triangulation, whose natural area is its Voronoi cell. + Cells which are unbounded (the mesh's convex hull) have no polygon and are omitted; bounded + cells are clipped to the bounding box of the mesh grid so a large boundary cell cannot stretch + the figure. + + Parameters + ---------- + mapper + A `Mapper` whose interpolator is Delaunay or K nearest neighbour. + pix_indexes + The indexes of the mesh pixels to return cells for. + + Returns + ------- + A list of ``(N, 2)`` arrays of ``(y, x)`` scaled coordinates, each a closed loop. + """ + voronoi = mapper.mesh_geometry.voronoi + + mesh_grid = np.asarray( + mapper.source_plane_mesh_grid.array + if hasattr(mapper.source_plane_mesh_grid, "array") + else mapper.source_plane_mesh_grid + ) + + y_min, y_max = float(np.min(mesh_grid[:, 0])), float(np.max(mesh_grid[:, 0])) + x_min, x_max = float(np.min(mesh_grid[:, 1])), float(np.max(mesh_grid[:, 1])) + + contours = [] + + for pix_index in pix_indexes: + region = voronoi.regions[voronoi.point_region[int(pix_index)]] + + if len(region) == 0 or -1 in region: + continue + + # `mesh_grid_xy` feeds scipy the mesh grid in its native (y, x) column order, so the + # Voronoi vertices come back in (y, x) too. + vertices = np.asarray(voronoi.vertices)[np.asarray(region, dtype=int)] + + vertices = np.stack( + [ + np.clip(vertices[:, 0], y_min, y_max), + np.clip(vertices[:, 1], x_min, x_max), + ], + axis=-1, + ) + + contours.append(np.vstack([vertices, vertices[0]])) + + return contours + + +def source_contours_from(mapper, pix_indexes) -> List[np.ndarray]: + """ + The source-plane cell boundaries of a group of mesh pixels, in scaled coordinates. + + Rectangular meshes return the four corners of each cell; Delaunay and K nearest neighbour + meshes return each vertex's Voronoi cell. + + Parameters + ---------- + mapper + The `Mapper` whose mesh the pixels belong to. + pix_indexes + The indexes of the mesh pixels forming the source-plane region. + + Returns + ------- + A list of ``(N, 2)`` arrays of ``(y, x)`` scaled coordinates, one closed loop per mesh pixel. + """ + from autoarray.inversion.mesh.interpolator.rectangular import ( + InterpolatorRectangular, + ) + from autoarray.inversion.mesh.interpolator.rectangular_uniform import ( + InterpolatorRectangularUniform, + ) + from autoarray.inversion.mesh.interpolator.delaunay import InterpolatorDelaunay + from autoarray.inversion.mesh.interpolator.knn import InterpolatorKNearestNeighbor + + pix_indexes = np.asarray(pix_indexes, dtype=int).ravel() + + if isinstance( + mapper.interpolator, (InterpolatorRectangular, InterpolatorRectangularUniform) + ): + y_edges, x_edges = _rectangular_edges_from(mapper=mapper) + + total_x = mapper.mesh_geometry.shape[1] + + contours = [] + + for pix_index in pix_indexes: + row = int(pix_index) // total_x + col = int(pix_index) % total_x + + y0, y1 = float(y_edges[row]), float(y_edges[row + 1]) + x0, x1 = float(x_edges[col]), float(x_edges[col + 1]) + + contours.append( + np.asarray([[y0, x0], [y0, x1], [y1, x1], [y1, x0], [y0, x0]]) + ) + + return contours + + if isinstance( + mapper.interpolator, (InterpolatorDelaunay, InterpolatorKNearestNeighbor) + ): + return _voronoi_contours_from(mapper=mapper, pix_indexes=pix_indexes) + + raise NotImplementedError( + f"Source-plane contours are not implemented for the interpolator " + f"{type(mapper.interpolator).__name__}." + ) + + +def mappings_from( + mapper, + pix_indexes, + weight_threshold: float = 0.0, + min_pixels: int = 1, + peak_values: Optional[List[float]] = None, +) -> List["Mapping"]: + """ + The `Mapping` objects of one or more groups of mesh pixels. + + Parameters + ---------- + mapper + The `Mapper` whose mappings between the data and the mesh are used. + pix_indexes + The mesh pixel indexes, either a flat sequence (one group) or a nested sequence (one group + per entry). + weight_threshold + A data pixel is in an image-plane region if its summed mapping weight exceeds this value. + min_pixels + Connected image-plane regions with fewer than this many pixels are discarded. + peak_values + The peak reconstructed value of each group, stored on the returned `Mapping` objects. `None` + leaves every `peak_value` as `None` (the bare-`Mapper` case, which has no reconstruction). + + Returns + ------- + One `Mapping` per group of mesh pixels. + """ + groups = pix_index_groups_from(pix_indexes=pix_indexes) + + mesh_grid = np.asarray( + mapper.source_plane_mesh_grid.array + if hasattr(mapper.source_plane_mesh_grid, "array") + else mapper.source_plane_mesh_grid + ) + + mappings = [] + + for i, group in enumerate(groups): + centre = np.mean(mesh_grid[group], axis=0) + + mappings.append( + Mapping( + pix_indexes=group, + source_contours=source_contours_from(mapper=mapper, pix_indexes=group), + source_centre=(float(centre[0]), float(centre[1])), + image_regions=image_regions_from( + mapper=mapper, + pix_indexes=group, + weight_threshold=weight_threshold, + min_pixels=min_pixels, + ), + peak_value=None if peak_values is None else float(peak_values[i]), + ) + ) + + return mappings diff --git a/autoarray/inversion/plot/inversion_plots.py b/autoarray/inversion/plot/inversion_plots.py index 9a9e32a5b..3c7445b5a 100644 --- a/autoarray/inversion/plot/inversion_plots.py +++ b/autoarray/inversion/plot/inversion_plots.py @@ -2,13 +2,22 @@ import logging import numpy as np from pathlib import Path -from typing import Optional, Union +from typing import List, Optional, Union from autonerves import conf from autoarray.inversion.mappers.abstract import Mapper from autoarray.plot.array import plot_array -from autoarray.plot.utils import subplots, numpy_grid, numpy_lines, numpy_positions, subplot_save, hide_unused_axes, conf_subplot_figsize, tight_layout +from autoarray.plot.utils import ( + subplots, + numpy_grid, + numpy_lines, + numpy_positions, + subplot_save, + hide_unused_axes, + conf_subplot_figsize, + tight_layout, +) from autoarray.inversion.plot.mapper_plots import plot_mapper from autoarray.structures.arrays.uniform_2d import Array2D @@ -235,6 +244,43 @@ def _recon_array(): subplot_save(fig, output_path, f"{output_filename}_{mapper_index}", output_format) +def _conf_mappings(key: str, default, old_key: Optional[str] = None): + """ + Read a `subplot_mappings` setting from the ``visualize/general.yaml`` ``inversion`` section. + + A config which predates the mappings rewrite has none of the new keys, so `old_key` names the + key it did have (``total_mappings_pixels``) and is read as a fallback for one release. When + neither key is present the documented `default` is used. + + Parameters + ---------- + key + The config key to read. + default + The value used when neither `key` nor `old_key` is in the config. + old_key + The superseded key read when `key` is absent. + + Returns + ------- + The config value. + """ + inversion = conf.instance["visualize"]["general"]["inversion"] + + try: + return inversion[key] + except KeyError: + pass + + if old_key is not None: + try: + return inversion[old_key] + except KeyError: + pass + + return default + + def subplot_mappings( inversion, pixelization_index: int = 0, @@ -247,9 +293,21 @@ def subplot_mappings( lines=None, grid=None, positions=None, + threshold: Optional[float] = None, + min_pixels: Optional[int] = None, + total_clumps: Optional[int] = None, + pix_indexes: Optional[List] = None, + weight_threshold: float = 0.0, + region_alpha: float = 0.25, ): """ - 2×2 subplot showing data, model image, reconstruction and unzoomed reconstruction. + 2x2 subplot showing how the brightest clumps of the source reconstruction map to the image-plane. + + The panels are the data with all other linear objects subtracted, the reconstructed image, the + source reconstruction zoomed to its brightest region, and the source reconstruction unzoomed. + Every source clump is drawn as a filled polygon on the two source-plane panels and its + image-plane regions (the multiple images it maps to) are drawn in the same colour on the two + image-plane panels, each labelled with the same number in both planes. Parameters ---------- @@ -269,22 +327,49 @@ def subplot_mappings( Apply log10 normalisation. mesh_grid, lines, grid, positions Optional overlays. + threshold + A source pixel joins a clump if its reconstructed value exceeds this fraction of the + reconstruction's maximum. ``None`` reads ``mappings_threshold`` from the visualize config. + min_pixels + Connected groups of source pixels smaller than this are not drawn. ``None`` reads + ``mappings_min_pixels`` from the visualize config. + total_clumps + The maximum number of clumps drawn, keeping the brightest. ``None`` reads + ``total_mappings`` from the visualize config. + pix_indexes + An explicit list of source pixel index groups which bypasses the clump finding, e.g. + ``[[0, 1], [10]]``. + weight_threshold + A data pixel is in an image-plane region if its summed mapping weight to the clump exceeds + this value. + region_alpha + The alpha of each drawn region's fill. """ mapper = inversion.cls_list_from(cls=Mapper)[pixelization_index] - try: - total_pixels = conf.instance["visualize"]["general"]["inversion"][ - "total_mappings_pixels" - ] - except Exception: - total_pixels = 10 + if threshold is None: + threshold = float(_conf_mappings(key="mappings_threshold", default=0.5)) + if min_pixels is None: + min_pixels = int(_conf_mappings(key="mappings_min_pixels", default=3)) + if total_clumps is None: + total_clumps = int( + _conf_mappings( + key="total_mappings", default=5, old_key="total_mappings_pixels" + ) + ) - pix_indexes = inversion.max_pixel_list_from( - total_pixels=total_pixels, - filter_neighbors=True, + mappings = inversion.mappings_from( mapper_index=pixelization_index, + weight_threshold=weight_threshold, + threshold=threshold, + min_pixels=min_pixels, + total_clumps=total_clumps, + pix_indexes=pix_indexes, ) - mapper.slim_indexes_for_pix_indexes(pix_indexes=pix_indexes) + + image_regions = [mapping.image_contours for mapping in mappings] + source_regions = [mapping.source_contours for mapping in mappings] + region_labels = [str(i + 1) for i in range(len(mappings))] fig, axes = subplots(2, 2, figsize=conf_subplot_figsize(2, 2)) axes = axes.flatten() @@ -305,6 +390,9 @@ def subplot_mappings( grid=grid, positions=positions, lines=lines, + regions=image_regions, + region_alpha=region_alpha, + region_labels=region_labels, ) except (AttributeError, KeyError): pass @@ -325,11 +413,15 @@ def subplot_mappings( grid=grid, positions=positions, lines=lines, + regions=image_regions, + region_alpha=region_alpha, + region_labels=region_labels, ) except (AttributeError, KeyError): pass pixel_values = inversion.reconstruction_dict[mapper] + plot_mapper( mapper, solution_vector=pixel_values, @@ -340,6 +432,9 @@ def subplot_mappings( zoom_to_brightest=True, mesh_grid=mesh_grid, lines=lines, + regions=source_regions, + region_alpha=region_alpha, + region_labels=region_labels, ) plot_mapper( mapper, @@ -351,6 +446,9 @@ def subplot_mappings( zoom_to_brightest=False, mesh_grid=mesh_grid, lines=lines, + regions=source_regions, + region_alpha=region_alpha, + region_labels=region_labels, ) hide_unused_axes(axes) @@ -402,9 +500,13 @@ def save_reconstruction_csv( ) noise_map = None - with open(output_path / f"source_plane_reconstruction_{i}.csv", mode="w", newline="") as f: + with open( + output_path / f"source_plane_reconstruction_{i}.csv", mode="w", newline="" + ) as f: writer = csv.writer(f) writer.writerow(["y", "x", "reconstruction", "noise_map"]) for j in range(len(x)): noise_value = float("nan") if noise_map is None else float(noise_map[j]) - writer.writerow([float(y[j]), float(x[j]), float(reconstruction[j]), noise_value]) + writer.writerow( + [float(y[j]), float(x[j]), float(reconstruction[j]), noise_value] + ) diff --git a/autoarray/inversion/plot/mapper_plots.py b/autoarray/inversion/plot/mapper_plots.py index 600851efa..c3028bedf 100644 --- a/autoarray/inversion/plot/mapper_plots.py +++ b/autoarray/inversion/plot/mapper_plots.py @@ -4,7 +4,14 @@ from autoarray.plot.array import plot_array from autoarray.plot.inversion import plot_inversion_reconstruction -from autoarray.plot.utils import subplots, numpy_grid, numpy_lines, subplot_save, conf_subplot_figsize, tight_layout +from autoarray.plot.utils import ( + subplots, + numpy_grid, + numpy_lines, + subplot_save, + conf_subplot_figsize, + tight_layout, +) logger = logging.getLogger(__name__) @@ -22,6 +29,10 @@ def plot_mapper( mesh_grid=None, lines=None, line_colors=None, + regions=None, + region_colors=None, + region_alpha: float = 0.25, + region_labels=None, title: str = "Pixelization Mesh (Source-Plane)", zoom_to_brightest: bool = True, zoom_extent_scale: float = 1.0, @@ -50,6 +61,15 @@ def plot_mapper( Mesh grid to overlay as scatter points. lines Lines to overlay. + regions + Source-plane regions to overlay as filled polygons, each a list of ``(N, 2)`` arrays of + ``(y, x)`` coordinates (e.g. the ``source_contours`` of a ``Mapping``). + region_colors + The colour of each region, cycled if there are more regions than colours. + region_alpha + The alpha of each region's fill. + region_labels + A text label drawn at the centre of each region. title Figure title. zoom_to_brightest @@ -71,6 +91,10 @@ def plot_mapper( zoom_extent_scale=zoom_extent_scale, lines=numpy_lines(lines), line_colors=line_colors, + regions=regions, + region_colors=region_colors, + region_alpha=region_alpha, + region_labels=region_labels, grid=numpy_grid(mesh_grid), output_path=output_path, output_filename=output_filename, @@ -90,6 +114,9 @@ def subplot_image_and_mapper( use_log10: bool = False, mesh_grid=None, lines=None, + regions=None, + region_colors=None, + region_labels=None, ): """ 1×2 subplot: image-plane image (left) and pixelization mesh (right). @@ -114,9 +141,24 @@ def subplot_image_and_mapper( Mesh grid to overlay on the reconstruction panel. lines Lines to overlay on both panels. + regions + A list of ``Mapping`` objects (e.g. from ``mapper.mappings_from(pix_indexes=...)``), drawn + as filled polygons in matched colours: each mapping's image-plane regions on the left panel + and its source-plane mesh pixels on the right panel. + region_colors + The colour of each region, cycled if there are more regions than colours. + region_labels + A text label drawn at the centre of each region, e.g. ``["1", "2"]``. """ fig, axes = subplots(1, 2, figsize=conf_subplot_figsize(1, 2)) + image_regions = ( + None if regions is None else [mapping.image_contours for mapping in regions] + ) + source_regions = ( + None if regions is None else [mapping.source_contours for mapping in regions] + ) + plot_array( image, ax=axes[0], @@ -124,6 +166,9 @@ def subplot_image_and_mapper( colormap=colormap, use_log10=use_log10, lines=lines, + regions=image_regions, + region_colors=region_colors, + region_labels=region_labels, ) plot_mapper( mapper, @@ -131,6 +176,9 @@ def subplot_image_and_mapper( use_log10=use_log10, mesh_grid=mesh_grid, lines=lines, + regions=source_regions, + region_colors=region_colors, + region_labels=region_labels, ax=axes[1], ) diff --git a/autoarray/plot/array.py b/autoarray/plot/array.py index b1de3efe2..b15466000 100644 --- a/autoarray/plot/array.py +++ b/autoarray/plot/array.py @@ -20,6 +20,7 @@ norm_from, _apply_contours, _conf_imshow_origin, + plot_regions, ) _zoom_array_2d = zoom_array @@ -44,6 +45,10 @@ def plot_array( patches: Optional[List] = None, fill_region: Optional[List] = None, contours: Optional[int] = None, + regions: Optional[List] = None, + region_colors: Optional[List] = None, + region_alpha: float = 0.25, + region_labels: Optional[List[str]] = None, # --- cosmetics -------------------------------------------------------------- title: str = "", xlabel: str = "", @@ -96,6 +101,18 @@ def plot_array( List of matplotlib ``Patch`` objects to draw over the image. fill_region List of two arrays ``[y1_arr, y2_arr]`` passed to ``ax.fill_between``. + regions + List of regions to overlay as filled polygons, each region being a list of ``(N, 2)`` + arrays of ``(y, x)`` coordinates (e.g. the ``contours`` of an ``ImageRegion``). Every + region is drawn in its own colour, so that the same colour identifies the same source + structure across figures. + region_colors + The colour of each region, cycled if there are more regions than colours. ``None`` uses + ``["r", "g", "b", "m", "c", "y"]``. + 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"]``. title Figure title string. xlabel, ylabel @@ -139,6 +156,7 @@ def plot_array( if colormap is None: from autoarray.plot.utils import _default_colormap + colormap = _default_colormap() # convert overlay params (safe for None and already-numpy inputs) @@ -199,6 +217,7 @@ def plot_array( if not is_rgb: from autoarray.plot.utils import _apply_colorbar + _apply_colorbar(im, ax, cb_unit=cb_unit, is_subplot=not owns_figure) # --- overlays -------------------------------------------------------------- @@ -238,17 +257,28 @@ def plot_array( pos = np.asarray(pos).reshape(-1, 2) ax.scatter(pos[:, 1], pos[:, 0], s=20, c=colors[i % len(colors)], zorder=5) - if lines is not None: for i, line in enumerate(lines): if line is not None and len(line) > 0: line = np.asarray(line).reshape(-1, 2) - color = line_colors[i] if (line_colors is not None and i < len(line_colors)) else None + color = ( + line_colors[i] + if (line_colors is not None and i < len(line_colors)) + else None + ) kw = {"linewidth": 1} if color is not None: kw["color"] = color ax.plot(line[:, 1], line[:, 0], **kw) + plot_regions( + ax, + regions=regions, + region_colors=region_colors, + region_alpha=region_alpha, + region_labels=region_labels, + ) + if vector_yx is not None: ax.quiver( vector_yx[:, 1], @@ -271,13 +301,17 @@ def plot_array( # Contours: auto-enabled for log10 plots; explicit int enables linear contours. if use_log10 or (contours is not None and contours > 0): _apply_contours( - ax, array, extent, + ax, + array, + extent, use_log10=use_log10, n=contours if (contours is not None and contours > 0) else None, ) # --- labels / ticks -------------------------------------------------------- - apply_labels(ax, title=title, xlabel=xlabel, ylabel=ylabel, is_subplot=not owns_figure) + apply_labels( + ax, title=title, xlabel=xlabel, ylabel=ylabel, is_subplot=not owns_figure + ) if extent is not None: apply_extent(ax, extent) diff --git a/autoarray/plot/inversion.py b/autoarray/plot/inversion.py index 984f9aceb..a99f465b6 100644 --- a/autoarray/plot/inversion.py +++ b/autoarray/plot/inversion.py @@ -8,7 +8,16 @@ import numpy as np -from autoarray.plot.utils import subplots, apply_extent, apply_labels, conf_figsize, save_figure, norm_from, _conf_imshow_origin +from autoarray.plot.utils import ( + subplots, + apply_extent, + apply_labels, + conf_figsize, + save_figure, + norm_from, + _conf_imshow_origin, + plot_regions, +) def plot_inversion_reconstruction( @@ -29,6 +38,10 @@ def plot_inversion_reconstruction( lines: Optional[List[np.ndarray]] = None, line_colors: Optional[List] = None, grid: Optional[np.ndarray] = None, + regions: Optional[List] = None, + region_colors: Optional[List] = None, + region_alpha: float = 0.25, + region_labels: Optional[List[str]] = None, # --- figure control (used only when ax is None) ----------------------------- figsize: Optional[Tuple[int, int]] = None, output_path: Optional[str] = None, @@ -64,6 +77,18 @@ def plot_inversion_reconstruction( Line overlays (e.g. critical curves). grid Scatter overlay (e.g. data-plane grid). + regions + List of source-plane regions to overlay as filled polygons, each region being a list of + ``(N, 2)`` arrays of ``(y, x)`` coordinates (e.g. the ``source_contours`` of a ``Mapping``). + Every region is drawn in its own colour, matching the colour its image-plane counterpart is + drawn with by :func:`~autoarray.plot.array.plot_array`. + region_colors + The colour of each region, cycled if there are more regions than colours. ``None`` uses + ``["r", "g", "b", "m", "c", "y"]``. + 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"]``. figsize, output_path, output_filename, output_format Figure output controls. """ @@ -78,6 +103,7 @@ def plot_inversion_reconstruction( if colormap is None: from autoarray.plot.utils import _default_colormap + colormap = _default_colormap() owns_figure = ax is None @@ -101,29 +127,51 @@ def plot_inversion_reconstruction( if isinstance( mapper.interpolator, (InterpolatorRectangular, InterpolatorRectangularUniform) ): - _plot_rectangular(ax, pixel_values, mapper, norm, colormap, extent, is_subplot=is_subplot) + _plot_rectangular( + ax, pixel_values, mapper, norm, colormap, extent, is_subplot=is_subplot + ) elif isinstance( mapper.interpolator, (InterpolatorDelaunay, InterpolatorKNearestNeighbor) ): - _plot_delaunay(ax, pixel_values, mapper, norm, colormap, extent, is_subplot=is_subplot) + _plot_delaunay( + ax, pixel_values, mapper, norm, colormap, extent, is_subplot=is_subplot + ) # --- overlays -------------------------------------------------------------- if lines is not None: for i, line in enumerate(lines): if line is not None and len(line) > 0: line = np.asarray(line).reshape(-1, 2) - color = line_colors[i] if (line_colors is not None and i < len(line_colors)) else None + color = ( + line_colors[i] + if (line_colors is not None and i < len(line_colors)) + else None + ) kw = {"linewidth": 1} if color is not None: kw["color"] = color ax.plot(line[:, 1], line[:, 0], **kw) + plot_regions( + ax, + regions=regions, + region_colors=region_colors, + region_alpha=region_alpha, + region_labels=region_labels, + ) + if grid is not None: ax.scatter(grid[:, 1], grid[:, 0], s=1, c="w", alpha=0.5) apply_extent(ax, extent) - apply_labels(ax, title=title, xlabel="" if is_subplot else xlabel, ylabel="" if is_subplot else ylabel, is_subplot=is_subplot) + apply_labels( + ax, + title=title, + xlabel="" if is_subplot else xlabel, + ylabel="" if is_subplot else ylabel, + is_subplot=is_subplot, + ) if owns_figure: save_figure( @@ -134,7 +182,9 @@ def plot_inversion_reconstruction( ) -def _plot_rectangular(ax, pixel_values, mapper, norm, colormap, extent, is_subplot=False): +def _plot_rectangular( + ax, pixel_values, mapper, norm, colormap, extent, is_subplot=False +): """Render a rectangular pixelization reconstruction onto *ax*. Uses ``imshow`` for uniform rectangular grids @@ -268,4 +318,5 @@ def _plot_delaunay(ax, pixel_values, mapper, norm, colormap, extent, is_subplot= tc = ax.tripcolor(x, y, vals, cmap=colormap, norm=norm) from autoarray.plot.utils import _apply_colorbar + _apply_colorbar(tc, ax, is_subplot=is_subplot) diff --git a/autoarray/plot/utils.py b/autoarray/plot/utils.py index 4e0793e81..3f8186d04 100644 --- a/autoarray/plot/utils.py +++ b/autoarray/plot/utils.py @@ -106,7 +106,12 @@ def zoom_array(array): except Exception: zoom_around_mask = False - if zoom_around_mask and hasattr(array, "mask") and hasattr(array.mask, "is_all_false") and not array.mask.is_all_false: + if ( + zoom_around_mask + and hasattr(array, "mask") + and hasattr(array.mask, "is_all_false") + and not array.mask.is_all_false + ): from autoarray.mask.derive.zoom_2d import Zoom2D return Zoom2D(mask=array.mask).array_2d_from(array=array, buffer=1) @@ -552,8 +557,12 @@ def apply_labels( """ if is_subplot: title_fs = conf_mat_plot_fontsize("title_subplot", default=20) - xlabel_fs = conf_mat_plot_fontsize("xlabel_subplot", default=conf_mat_plot_fontsize("xlabel", default=14)) - ylabel_fs = conf_mat_plot_fontsize("ylabel_subplot", default=conf_mat_plot_fontsize("ylabel", default=14)) + xlabel_fs = conf_mat_plot_fontsize( + "xlabel_subplot", default=conf_mat_plot_fontsize("xlabel", default=14) + ) + ylabel_fs = conf_mat_plot_fontsize( + "ylabel_subplot", default=conf_mat_plot_fontsize("ylabel", default=14) + ) xticks_fs = conf_mat_plot_fontsize("xticks_subplot", default=18) yticks_fs = conf_mat_plot_fontsize("yticks_subplot", default=18) else: @@ -609,6 +618,7 @@ def save_figure( if dpi is None: from autonerves import conf + dpi = int(conf.instance["visualize"]["general"]["general"]["dpi"]) if _output_mode_save(fig, filename): @@ -635,9 +645,7 @@ def save_figure( pad_inches=0.1, ) except Exception as exc: - logger.warning( - f"save_figure: could not save {filename}.{fmt}: {exc}" - ) + logger.warning(f"save_figure: could not save {filename}.{fmt}: {exc}") plt.close(fig) @@ -672,6 +680,7 @@ def plot_visibilities_1d(vis, ax, title: str = "") -> None: def _conf_colorbar(key: str, default): try: from autonerves import conf + return conf.instance["visualize"]["general"]["colorbar"][key] except Exception: return default @@ -682,6 +691,7 @@ def _colorbar_tick_values(norm) -> Optional[List[float]]: if norm is None or norm.vmin is None or norm.vmax is None: return None import matplotlib.colors as mcolors + lo, hi = float(norm.vmin), float(norm.vmax) if isinstance(norm, mcolors.LogNorm): mid = 10 ** ((np.log10(lo) + np.log10(hi)) / 2.0) @@ -725,7 +735,9 @@ def _fmt_tick(v: float) -> str: return f"{v:.2f}" -def _colorbar_tick_labels(tick_values: List[float], cb_unit: Optional[str] = None) -> List[str]: +def _colorbar_tick_labels( + tick_values: List[float], cb_unit: Optional[str] = None +) -> List[str]: """Format tick values, appending *cb_unit* to the middle label. All three labels use a consistent notation style: if any tick is rendered @@ -740,6 +752,7 @@ def _colorbar_tick_labels(tick_values: List[float], cb_unit: Optional[str] = Non if cb_unit is None: try: from autonerves import conf + cb_unit = conf.instance["visualize"]["general"]["units"]["cb_unit"] except Exception: cb_unit = "" @@ -840,6 +853,7 @@ def _apply_contours( """ try: from autonerves import conf + _c = conf.instance["visualize"]["general"]["contour"] total = int(n if n is not None else _c.get("total_contours", 10)) include_values = bool(_c.get("include_values", True)) @@ -851,7 +865,10 @@ def _apply_contours( if use_log10: try: from autonerves import conf - log10_min = float(conf.instance["visualize"]["general"]["general"]["log10_min_value"]) + + log10_min = float( + conf.instance["visualize"]["general"]["general"]["log10_min_value"] + ) except Exception: log10_min = 1.0e-4 @@ -867,7 +884,9 @@ def _apply_contours( total, ) else: - levels = np.linspace(float(np.nanmin(array)), float(np.nanmax(array)), total) + levels = np.linspace( + float(np.nanmin(array)), float(np.nanmax(array)), total + ) # Build explicit coordinate grids so the contours align with imshow. # With origin="upper" row 0 maps to ymax (Y decreases across rows); @@ -1035,6 +1054,7 @@ def _conf_log10_min_value() -> float: """Return the log10 colour-scale floor from config (``log10_min_value``).""" try: from autonerves import conf + general = conf.instance["visualize"]["general"]["general"] return float(general["log10_min_value"]) except Exception: @@ -1044,6 +1064,7 @@ def _conf_log10_min_value() -> float: def _conf_ticks(key: str, default: float) -> float: try: from autonerves import conf + return float(conf.instance["visualize"]["general"]["ticks"][key]) except Exception: return default @@ -1052,6 +1073,7 @@ def _conf_ticks(key: str, default: float) -> float: def _conf_ticks_flag(key: str, default: bool) -> bool: try: from autonerves import conf + return bool(conf.instance["visualize"]["general"]["ticks"][key]) except Exception: return default @@ -1104,7 +1126,7 @@ def _fmt_minus(label: str) -> str: # render with a proper minus sign rather than a hyphen. return label.replace("-", "−") if minus_in_math else label - labels = [f'{v:g}' for v in ticks] + labels = [f"{v:g}" for v in ticks] if any("." in label for label in labels): labels = [ label if "." in label else f"{float(v):.1f}" @@ -1152,3 +1174,81 @@ def apply_extent( # drops below the tick line. Centre them on the tick. for label in ax.get_yticklabels(): label.set_verticalalignment("center") + + +REGION_COLORS_DEFAULT = ["r", "g", "b", "m", "c", "y"] + + +def plot_regions( + ax, + regions, + region_colors: Optional[List] = None, + region_alpha: float = 0.25, + region_labels: Optional[List[str]] = None, +) -> None: + """ + Draw filled and outlined polygon regions onto *ax*, one colour per region. + + Every region is drawn as one or more closed polygons in the plot's coordinate space (e.g. + arc-seconds), so the overlay stays aligned with the image however it is zoomed or cropped. + + Parameters + ---------- + ax + The matplotlib ``Axes`` to draw onto. + regions + A list of regions, where each region is a list of ``(N, 2)`` arrays of ``(y, x)`` + coordinates (a single ``(N, 2)`` array is also accepted for a one-polygon region). + region_colors + The colour of each region, cycled if there are more regions than colours. ``None`` uses the + default cycle ``["r", "g", "b", "m", "c", "y"]``. + 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. + """ + if regions is None: + return + + colors = region_colors if region_colors is not None else REGION_COLORS_DEFAULT + + for i, region in enumerate(regions): + if region is None: + continue + + if isinstance(region, np.ndarray) and region.ndim == 2: + polygons = [region] + else: + polygons = list(region) + + color = colors[i % len(colors)] + + points = [] + + for polygon in polygons: + polygon = np.asarray(polygon, dtype=float).reshape(-1, 2) + + if polygon.shape[0] < 3: + continue + + points.append(polygon) + + ax.fill( + polygon[:, 1], polygon[:, 0], color=color, alpha=region_alpha, zorder=4 + ) + 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, + ) diff --git a/test_autoarray/inversion/inversion/test_abstract.py b/test_autoarray/inversion/inversion/test_abstract.py index e5c8f34de..3b10a6588 100644 --- a/test_autoarray/inversion/inversion/test_abstract.py +++ b/test_autoarray/inversion/inversion/test_abstract.py @@ -1087,3 +1087,110 @@ def test__reconstruction_covariance_matrix__non_finite_entry_only_raises_when_th with pytest.raises(np.linalg.LinAlgError, match="non-finite"): _zeroed_pixel_inversion(kept_nan).reconstruction_covariance_matrix + + +def _inversion_with_reconstruction(monkeypatch, inversion, reconstruction): + """Force an inversion's reconstruction to a chosen vector, so clump finding can be tested.""" + mapper = inversion.cls_list_from(cls=aa.Mapper)[0] + + monkeypatch.setattr( + type(inversion), + "reconstruction_dict", + property(lambda self: {mapper: np.asarray(reconstruction)}), + ) + + return inversion + + +def test__source_clumps_from__two_peaks_give_two_clumps( + rectangular_inversion_7x7_3x3, monkeypatch +): + # A 3x3 mesh with bright opposite corners, whose bright pixels do not touch. + inversion = _inversion_with_reconstruction( + monkeypatch, + rectangular_inversion_7x7_3x3, + [1.0, 0.6, 0.1, 0.6, 0.1, 0.1, 0.1, 0.6, 1.0], + ) + + clumps = inversion.source_clumps_from(threshold=0.5, min_pixels=1) + + assert len(clumps) == 2 + assert set(clumps[0].tolist()) == {0, 1, 3} + assert set(clumps[1].tolist()) == {7, 8} + + +def test__source_clumps_from__low_threshold_merges_the_peaks_into_one_clump( + rectangular_inversion_7x7_3x3, monkeypatch +): + inversion = _inversion_with_reconstruction( + monkeypatch, + rectangular_inversion_7x7_3x3, + [1.0, 0.6, 0.1, 0.6, 0.1, 0.1, 0.1, 0.6, 1.0], + ) + + clumps = inversion.source_clumps_from(threshold=0.05, min_pixels=1) + + assert len(clumps) == 1 + assert clumps[0].shape[0] == 9 + + +def test__source_clumps_from__min_pixels_and_total_clumps_filter_the_clumps( + rectangular_inversion_7x7_3x3, monkeypatch +): + inversion = _inversion_with_reconstruction( + monkeypatch, + rectangular_inversion_7x7_3x3, + [1.0, 0.6, 0.1, 0.6, 0.1, 0.1, 0.1, 0.6, 1.0], + ) + + assert len(inversion.source_clumps_from(threshold=0.5, min_pixels=3)) == 1 + assert len(inversion.source_clumps_from(threshold=0.5, min_pixels=4)) == 0 + assert ( + len(inversion.source_clumps_from(threshold=0.5, min_pixels=1, total_clumps=1)) + == 1 + ) + + +def test__source_clumps_from__pix_indexes_bypasses_the_clump_finding( + rectangular_inversion_7x7_3x3, +): + clumps = rectangular_inversion_7x7_3x3.source_clumps_from(pix_indexes=[[0, 1], [8]]) + + assert len(clumps) == 2 + assert (clumps[0] == np.array([0, 1])).all() + assert (clumps[1] == np.array([8])).all() + + +def test__source_clumps_from__non_positive_reconstruction_gives_no_clumps( + rectangular_inversion_7x7_3x3, monkeypatch +): + inversion = _inversion_with_reconstruction( + monkeypatch, rectangular_inversion_7x7_3x3, -1.0 * np.ones(9) + ) + + assert inversion.source_clumps_from(min_pixels=1) == [] + + +def test__mappings_from__pairs_each_clump_with_its_image_regions( + rectangular_inversion_7x7_3x3, monkeypatch +): + inversion = _inversion_with_reconstruction( + monkeypatch, + rectangular_inversion_7x7_3x3, + [1.5, 0.6, 0.1, 0.6, 0.1, 0.1, 0.1, 1.2, 2.0], + ) + + mappings = inversion.mappings_from(threshold=0.5, min_pixels=1) + + assert len(mappings) == 2 + + # Ordered brightest first, so the clump containing mesh pixel 8 (value 2.0) comes first. + assert mappings[0].peak_value == pytest.approx(2.0) + assert set(mappings[0].pix_indexes.tolist()) == {7, 8} + assert mappings[1].peak_value == pytest.approx(1.5) + assert set(mappings[1].pix_indexes.tolist()) == {0} + + for mapping in mappings: + assert len(mapping.source_contours) == mapping.pix_indexes.shape[0] + assert len(mapping.image_regions) > 0 + assert len(mapping.source_centre) == 2 diff --git a/test_autoarray/inversion/mappings/__init__.py b/test_autoarray/inversion/mappings/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/test_autoarray/inversion/mappings/test_mapping.py b/test_autoarray/inversion/mappings/test_mapping.py new file mode 100644 index 000000000..7bcb2cdca --- /dev/null +++ b/test_autoarray/inversion/mappings/test_mapping.py @@ -0,0 +1,237 @@ +import numpy as np +import pytest + +import autoarray as aa +from autoarray.inversion.linear_obj.neighbors import Neighbors +from autoarray.inversion.mappings.mapping import ( + connected_components_from, + contours_from_bool_native, + image_regions_from, + image_regions_from_slim_mask, + pix_index_groups_from, + source_contours_from, +) + + +class MockMapper: + """A mapper stripped down to what the image-plane region helpers read from it.""" + + def __init__(self, mapping_matrix, mask): + self.mapping_matrix = mapping_matrix + self.mask = mask + + +def make_mask_5x5(): + return aa.Mask2D.all_false(shape_native=(5, 5), pixel_scales=1.0) + + +def make_bimodal_mapper(): + """ + A mapper whose single mesh pixel maps to two disconnected 2x2 blocks of a 5x5 mask -- the + two "multiple images" of one source region. + """ + mask = make_mask_5x5() + + mapping_matrix = np.zeros(shape=(25, 2)) + mapping_matrix[[0, 1, 5, 6], 0] = 1.0 + mapping_matrix[[18, 19, 23, 24], 0] = 1.0 + mapping_matrix[[12], 1] = 1.0 + + return MockMapper(mapping_matrix=mapping_matrix, mask=mask) + + +def test__connected_components_from__splits_a_disconnected_index_set(): + # A chain 0 - 1 - 2, plus an isolated pixel 3. + neighbors = Neighbors( + arr=np.array([[1, -1], [0, 2], [1, -1], [-1, -1]]), + sizes=np.array([1, 2, 1, 0]), + ) + + components = connected_components_from(indexes=[0, 1, 2, 3], neighbors=neighbors) + + assert len(components) == 2 + assert (components[0] == np.array([0, 1, 2])).all() + assert (components[1] == np.array([3])).all() + + +def test__connected_components_from__only_walks_through_indexes_in_the_input(): + neighbors = Neighbors( + arr=np.array([[1, -1], [0, 2], [1, -1], [-1, -1]]), + sizes=np.array([1, 2, 1, 0]), + ) + + # Pixel 1 is not in the input, so 0 and 2 are not connected through it. + components = connected_components_from(indexes=[0, 2], neighbors=neighbors) + + assert len(components) == 2 + + +def test__connected_components_from__ignores_entries_beyond_the_neighbor_size(): + # The rectangular mesh leaves stale (non -1) values beyond `sizes`, which must be ignored. + neighbors = Neighbors( + arr=np.array([[1, 3], [0, 3], [3, 3]]), + sizes=np.array([1, 1, 0]), + ) + + components = connected_components_from(indexes=[0, 1, 2], neighbors=neighbors) + + assert len(components) == 2 + assert (components[0] == np.array([0, 1])).all() + assert (components[1] == np.array([2])).all() + + +def test__pix_index_groups_from__flat_and_nested_inputs(): + assert len(pix_index_groups_from(pix_indexes=[0, 1, 2])) == 1 + assert len(pix_index_groups_from(pix_indexes=[[0, 1], [2]])) == 2 + assert pix_index_groups_from(pix_indexes=[]) == [] + + +def test__contours_from_bool_native__single_block_is_one_closed_loop(): + bool_native = np.zeros(shape=(5, 5), dtype=bool) + bool_native[0:2, 0:2] = True + + contours = contours_from_bool_native( + bool_native=bool_native, geometry=make_mask_5x5().geometry + ) + + assert len(contours) == 1 + + contour = contours[0] + + # One point per pixel edge of the block's boundary (8 edges), closed back onto its start. + assert contour.shape == (9, 2) + assert contour[0] == pytest.approx(contour[-1]) + + # The block's corners in scaled units: y from +2.5 down to +0.5, x from -2.5 to -0.5. + assert contour[:, 0].max() == pytest.approx(2.5) + assert contour[:, 0].min() == pytest.approx(0.5) + assert contour[:, 1].max() == pytest.approx(-0.5) + assert contour[:, 1].min() == pytest.approx(-2.5) + + +def test__image_regions_from__bimodal_mapping_matrix_gives_two_regions(): + mapper = make_bimodal_mapper() + + regions = image_regions_from(mapper=mapper, pix_indexes=[0], min_pixels=1) + + assert len(regions) == 2 + + assert (regions[0].slim_indexes == np.array([0, 1, 5, 6])).all() + assert (regions[1].slim_indexes == np.array([18, 19, 23, 24])).all() + + for region in regions: + assert len(region.contours) == 1 + assert region.contours[0][0] == pytest.approx(region.contours[0][-1]) + assert region.area() == pytest.approx(4.0) + + +def test__image_regions_from__min_pixels_discards_small_regions(): + mapper = make_bimodal_mapper() + + # Mesh pixel 1 maps to a single data pixel only. + assert len(image_regions_from(mapper=mapper, pix_indexes=[1], min_pixels=1)) == 1 + assert len(image_regions_from(mapper=mapper, pix_indexes=[1], min_pixels=2)) == 0 + + +def test__image_regions_from__weight_threshold_removes_weakly_mapped_pixels(): + mask = make_mask_5x5() + + mapping_matrix = np.zeros(shape=(25, 1)) + mapping_matrix[[0, 1, 5, 6], 0] = 1.0 + mapping_matrix[[24], 0] = 0.1 + + mapper = MockMapper(mapping_matrix=mapping_matrix, mask=mask) + + assert len(image_regions_from(mapper=mapper, pix_indexes=[0], min_pixels=1)) == 2 + assert ( + len( + image_regions_from( + mapper=mapper, pix_indexes=[0], weight_threshold=0.5, min_pixels=1 + ) + ) + == 1 + ) + + +def test__image_regions_from_slim_mask__region_centre_and_flux_methods(): + mask = make_mask_5x5() + + slim_bool = np.zeros(25, dtype=bool) + slim_bool[[0, 1, 5, 6]] = True + + regions = image_regions_from_slim_mask(mask=mask, slim_bool=slim_bool, min_pixels=1) + + assert len(regions) == 1 + + region = regions[0] + + assert region.centre == pytest.approx((1.5, -1.5)) + + values = np.zeros(shape=(5, 5)) + values[0, 1] = 4.0 + values[1, 0] = 1.0 + + array = aa.Array2D(values=values, mask=mask) + + assert region.flux_from(array=array) == pytest.approx(5.0) + assert region.brightest_coordinate_from(array=array) == pytest.approx((2.0, -1.0)) + # Flux-weighted: (4 * (2.0, -1.0) + 1 * (1.0, -2.0)) / 5. + assert region.centroid_from(array=array) == pytest.approx((1.8, -1.2)) + + +def test__source_contours_from__rectangular_returns_one_closed_cell_per_pixel( + rectangular_mapper_7x7_3x3, +): + contours = source_contours_from( + mapper=rectangular_mapper_7x7_3x3, pix_indexes=[0, 4] + ) + + assert len(contours) == 2 + + for contour in contours: + assert contour.shape == (5, 2) + assert contour[0] == pytest.approx(contour[-1]) + + # The central mesh pixel's cell is centred on the mesh origin. + assert np.mean(contours[1][:-1, 0]) == pytest.approx(0.0, abs=1.0e-6) + assert np.mean(contours[1][:-1, 1]) == pytest.approx(0.0, abs=1.0e-6) + + +def test__source_contours_from__delaunay_returns_closed_voronoi_cells( + delaunay_mapper_9_3x3, +): + contours = source_contours_from( + mapper=delaunay_mapper_9_3x3, pix_indexes=list(range(9)) + ) + + # Unbounded cells on the mesh's convex hull have no polygon, so fewer cells than pixels. + assert 0 < len(contours) <= 9 + + for contour in contours: + assert contour.shape[1] == 2 + assert contour[0] == pytest.approx(contour[-1]) + + +def test__mapper_mappings_from__pairs_source_pixels_with_image_regions( + rectangular_mapper_7x7_3x3, +): + mappings = rectangular_mapper_7x7_3x3.mappings_from(pix_indexes=[[0, 1], [8]]) + + assert len(mappings) == 2 + + assert (mappings[0].pix_indexes == np.array([0, 1])).all() + assert mappings[0].peak_value is None + assert len(mappings[0].source_contours) == 2 + assert len(mappings[0].image_regions) > 0 + assert len(mappings[0].image_contours) == sum( + len(region.contours) for region in mappings[0].image_regions + ) + + +def test__mapper_mappings_from__flat_pix_indexes_is_one_mapping( + rectangular_mapper_7x7_3x3, +): + mappings = rectangular_mapper_7x7_3x3.mappings_from(pix_indexes=[0, 1]) + + assert len(mappings) == 1 + assert (mappings[0].pix_indexes == np.array([0, 1])).all() diff --git a/test_autoarray/inversion/plot/test_inversion_plotters.py b/test_autoarray/inversion/plot/test_inversion_plotters.py index ccb0aa595..c43b4e128 100644 --- a/test_autoarray/inversion/plot/test_inversion_plotters.py +++ b/test_autoarray/inversion/plot/test_inversion_plotters.py @@ -137,3 +137,37 @@ def test__save_reconstruction_csv__singular_curvature_reg_matrix( assert np.isfinite(float(row[0])) assert np.isfinite(float(row[1])) assert np.isnan(float(row[3])) + + +def test__subplot_mappings__delaunay_inversion( + delaunay_inversion_9_3x3, + plot_path, + plot_patch, +): + aplt.subplot_mappings( + inversion=delaunay_inversion_9_3x3, + pixelization_index=0, + min_pixels=1, + output_path=plot_path, + output_filename="mappings_delaunay", + output_format="png", + ) + + assert str(Path(plot_path) / "mappings_delaunay_0.png") in plot_patch.paths + + +def test__subplot_mappings__pix_indexes_bypasses_the_clump_finding( + rectangular_inversion_7x7_3x3, + plot_path, + plot_patch, +): + aplt.subplot_mappings( + inversion=rectangular_inversion_7x7_3x3, + pixelization_index=0, + pix_indexes=[[0, 1], [8]], + output_path=plot_path, + output_filename="mappings_pix_indexes", + output_format="png", + ) + + assert str(Path(plot_path) / "mappings_pix_indexes_0.png") in plot_patch.paths diff --git a/test_autoarray/inversion/plot/test_mapper_plotters.py b/test_autoarray/inversion/plot/test_mapper_plotters.py index 3e14ccb99..72e3bf00a 100644 --- a/test_autoarray/inversion/plot/test_mapper_plotters.py +++ b/test_autoarray/inversion/plot/test_mapper_plotters.py @@ -58,3 +58,24 @@ def test__subplot_image_and_mapper( output_format="png", ) assert str(Path(plot_path) / "image_and_mapper.png") in plot_patch.paths + + +def test__subplot_image_and_mapper__with_regions( + imaging_7x7, + rectangular_mapper_7x7_3x3, + plot_path, + plot_patch, +): + mappings = rectangular_mapper_7x7_3x3.mappings_from(pix_indexes=[[0, 1], [8]]) + + aplt.subplot_image_and_mapper( + mapper=rectangular_mapper_7x7_3x3, + image=imaging_7x7.data, + regions=mappings, + region_labels=["1", "2"], + output_path=plot_path, + output_filename="image_and_mapper_regions", + output_format="png", + ) + + assert str(Path(plot_path) / "image_and_mapper_regions.png") in plot_patch.paths diff --git a/test_autoarray/plot/test_utils.py b/test_autoarray/plot/test_utils.py index 30399267a..5d4334cda 100644 --- a/test_autoarray/plot/test_utils.py +++ b/test_autoarray/plot/test_utils.py @@ -74,7 +74,7 @@ def test_arcsec_labels_minus_in_math(): # Suffix format: the ASCII hyphen (U+002D) becomes the math minus (U+2212). ticks["symbol_over_decimal"] = False - assert _arcsec_labels([-2.1, -0.044, 2.0]) == ["−2.1\"", "−0.044\"", "2.0\""] + assert _arcsec_labels([-2.1, -0.044, 2.0]) == ['−2.1"', '−0.044"', '2.0"'] # Symbol-over-decimal format also uses the math minus. ticks["symbol_over_decimal"] = True @@ -471,3 +471,37 @@ def test__invalid_value_raises(self): general["imshow_origin"] = original assert "sideways" in str(exc_info.value) + + +def test__plot_array__regions_overlay_is_output(array_2d_7x7, plot_path, plot_patch): + """A `regions` polygon is arcsec-space, so it survives the `zoom_array` crop `plot_array` applies.""" + region = np.array([[1.0, -1.0], [1.0, 1.0], [-1.0, 1.0], [-1.0, -1.0], [1.0, -1.0]]) + + aplt.plot_array( + array=array_2d_7x7, + regions=[[region]], + region_labels=["1"], + output_path=plot_path, + output_filename="array_regions", + output_format="png", + ) + + assert str(Path(plot_path) / "array_regions.png") in plot_patch.paths + + +def test__plot_inversion_reconstruction__regions_overlay_is_output( + rectangular_mapper_7x7_3x3, plot_path, plot_patch +): + mappings = rectangular_mapper_7x7_3x3.mappings_from(pix_indexes=[[0, 1], [8]]) + + aplt.plot_inversion_reconstruction( + pixel_values=np.ones(9), + mapper=rectangular_mapper_7x7_3x3, + regions=[mapping.source_contours for mapping in mappings], + region_labels=["1", "2"], + output_path=plot_path, + output_filename="reconstruction_regions", + output_format="png", + ) + + assert str(Path(plot_path) / "reconstruction_regions.png") in plot_patch.paths