From cef287163f70e8becff045d0b3e22caac92796bb Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 15 May 2026 15:55:45 +0100 Subject: [PATCH] fix(hilbert): support offset-centre circular masks Previously Hilbert image-mesh raised PixelizationException for offset-centre circular masks whose centres didn't land on a pixel half-integer, and would have silently sampled the wrong region of the image for those that slipped through. This commit fixes both halves of the bug. - Mask2D.is_circular rewritten to be quantization-robust: bbox square within one pixel, geometric centre pixel unmasked, and a reference circular mask reconstructed from the inferred centre+radius must match the input within rim quantization tolerance. The new predicate correctly accepts offset circles and rejects annular and square masks that previously false-passed. - Mask2D.circular_radius now derives the radius from the bbox extent rather than the central row count. Numerically identical for centred masks; correct for offset centres. - hilbert.image_and_grid_from now translates Hilbert sample points by mask_centre, so the curve actually probes the masked region of the image. For masks centred at the origin (every current test/smoke case) the shift is a no-op and likelihoods are bit-identical. Co-Authored-By: Claude Opus 4.7 (1M context) --- .../inversion/mesh/image_mesh/hilbert.py | 10 ++- autoarray/mask/mask_2d.py | 81 ++++++++++++++----- .../pixelization/image_mesh/test_hilbert.py | 28 +++++++ test_autoarray/mask/test_mask_2d.py | 60 ++++++++++++++ 4 files changed, 160 insertions(+), 19 deletions(-) diff --git a/autoarray/inversion/mesh/image_mesh/hilbert.py b/autoarray/inversion/mesh/image_mesh/hilbert.py index 945bd860b..85f1a2295 100644 --- a/autoarray/inversion/mesh/image_mesh/hilbert.py +++ b/autoarray/inversion/mesh/image_mesh/hilbert.py @@ -118,10 +118,16 @@ def grid_hilbert_order_from(length, mask_radius): return x1d_hb, y1d_hb -def image_and_grid_from(image, mask, mask_radius, pixel_scales, hilbert_length): +def image_and_grid_from( + image, mask, mask_radius, pixel_scales, hilbert_length, mask_centre=(0.0, 0.0) +): """ This code will create a grid in Hilbert space-filling curve order and an interpolated hyper image associated to that grid. + + The Hilbert curve is generated around the origin and filtered to a disc of radius ``mask_radius``, then translated + by ``mask_centre`` so the sampling lands inside the mask's circle for offset-centre masks. For masks centred on + the origin the translation is a no-op and behaviour is unchanged. """ from scipy.interpolate import griddata @@ -145,6 +151,7 @@ def image_and_grid_from(image, mask, mask_radius, pixel_scales, hilbert_length): grid_hb = np.stack((y1d_hb, x1d_hb), axis=-1) grid_hb_radius = np.sqrt(grid_hb[:, 0] ** 2.0 + grid_hb[:, 1] ** 2.0) new_grid = grid_hb[grid_hb_radius <= mask_radius] + new_grid = new_grid + np.asarray(mask_centre) new_img = griddata( points=grid, @@ -271,6 +278,7 @@ def image_plane_mesh_grid_from( mask_radius=mask.circular_radius, pixel_scales=mask.pixel_scales, hilbert_length=193, + mask_centre=mask.mask_centre, ) weight_map = self.weight_map_from(adapt_data=adapt_data_hb) diff --git a/autoarray/mask/mask_2d.py b/autoarray/mask/mask_2d.py index 3fed1b2f1..66c087498 100644 --- a/autoarray/mask/mask_2d.py +++ b/autoarray/mask/mask_2d.py @@ -938,10 +938,23 @@ def resized_from(self, new_shape, pad_value: int = 0.0) -> Mask2D: @property def is_circular(self) -> bool: """ - Returns whether the mask is circular or not. + Returns whether the unmasked region of the mask is a filled circular disc. - This is performed by taking the central row and column of the mask (based on the mask centre) and counting - the number of unmasked pixels. If the number of unmasked pixels is the same, the mask is circular. + The check is robust to circles whose centre is offset from the coordinate origin, including offsets that + do not align with pixel boundaries. It rejects annular, square and other non-disc shapes. + + The algorithm: + + 1) The bounding box of unmasked pixels must be square to within one pixel. Offset centres that fall between + pixel centres can produce a one-pixel asymmetry in the bounding box, which is allowed; anything larger + indicates a non-circular shape such as an ellipse. Tiny masks (bounding box <= 2 pixels) require an exactly + square bounding box, since the one-pixel slack would otherwise admit 1x2 strips. + 2) The pixel at the geometric centre of the bounding box must be unmasked. This rejects annular masks whose + inner hole is at least one pixel wide. + 3) The centre and radius of the unmasked region are inferred from the bounding box and pixel count, and a + reference circular mask is built with those parameters. The input mask must match the reference within a + small number of pixels (tolerance scales with mask area to absorb rim quantization). This rejects squares, + crosses and tight annuli that slipped past the earlier checks. This function does not support rectangular masks and an exception will be raised if the pixel scales in each direction are different. @@ -955,22 +968,55 @@ def is_circular(self) -> bool: """ ) - pixel_coordinates_2d = self.geometry.pixel_coordinates_2d_from( - scaled_coordinates_2d=self.mask_centre + where = np.where(np.invert(self.array)) + if where[0].size == 0: + return False + + y_min, y_max = int(where[0].min()), int(where[0].max()) + x_min, x_max = int(where[1].min()), int(where[1].max()) + y_extent = y_max - y_min + 1 + x_extent = x_max - x_min + 1 + + if abs(y_extent - x_extent) > 1: + return False + if max(y_extent, x_extent) <= 2 and y_extent != x_extent: + return False + + cy = (y_max + y_min) // 2 + cx = (x_max + x_min) // 2 + if bool(self.array[cy, cx]): + return False + + actual_area = int(where[0].size) + inferred_radius_pix = np.sqrt(actual_area / np.pi) + inferred_radius = inferred_radius_pix * self.pixel_scales[0] + + y_centre_scaled = ( + 0.5 * (self.shape_native[0] - 1) - 0.5 * (y_min + y_max) + ) * self.pixel_scales[0] + x_centre_scaled = ( + 0.5 * (x_min + x_max) - 0.5 * (self.shape_native[1] - 1) + ) * self.pixel_scales[1] + + expected = mask_2d_util.mask_2d_circular_from( + shape_native=self.shape_native, + pixel_scales=self.pixel_scales, + radius=inferred_radius, + centre=(y_centre_scaled, x_centre_scaled), ) - central_row_pixels = sum(np.invert(self[pixel_coordinates_2d[0], :])) - central_column_pixels = sum(np.invert(self[:, pixel_coordinates_2d[1]])) - - return central_row_pixels == central_column_pixels + diff_count = int(np.sum(self.array != expected)) + tolerance = max(2, int(0.1 * actual_area)) + return diff_count <= tolerance @property def circular_radius(self) -> float: """ Returns the radius in scaled units of a circular mask. - This is performed by taking the central row of the mask (based on the mask centre) and counting the number of - unmasked pixels. The radius is then half the number of unmasked pixels times the pixel scale. + The radius is computed from the bounding box of the unmasked region, taking the larger of the y and x extents + as the diameter in pixels. This is robust to offset centres that fall between pixel boundaries (where the + bounding box can be asymmetric by one pixel). The mask is first checked that it is circular using the `is_circular` property, with an exception raised if it is not. @@ -987,15 +1033,14 @@ def circular_radius(self) -> float: raise exc.MaskException( """ A circular radius can only be computed for a circular mask. - + The `is_circular` property of this mask has returned False, indicating the mask is not circular. """ ) - pixel_coordinates_2d = self.geometry.pixel_coordinates_2d_from( - scaled_coordinates_2d=self.mask_centre - ) - - central_row_pixels = sum(np.invert(self[pixel_coordinates_2d[0], :])) + where = np.where(np.invert(self.array)) + y_extent = int(where[0].max() - where[0].min() + 1) + x_extent = int(where[1].max() - where[1].min() + 1) + diameter = max(y_extent, x_extent) - return central_row_pixels * self.pixel_scales[0] / 2.0 + return diameter * self.pixel_scales[0] / 2.0 diff --git a/test_autoarray/inversion/pixelization/image_mesh/test_hilbert.py b/test_autoarray/inversion/pixelization/image_mesh/test_hilbert.py index 5a85904bd..8654f23df 100644 --- a/test_autoarray/inversion/pixelization/image_mesh/test_hilbert.py +++ b/test_autoarray/inversion/pixelization/image_mesh/test_hilbert.py @@ -1,3 +1,4 @@ +import numpy as np import pytest import autoarray as aa @@ -22,3 +23,30 @@ def test__image_plane_mesh_grid_from(): [-1.02590674, -1.70984456], 1.0e-4, ) + + +def test__image_plane_mesh_grid_from__offset_centre__points_inside_mask_circle(): + centre = (0.5, 0.3) + radius = 0.5 + + mask = aa.Mask2D.circular( + shape_native=(40, 40), + radius=radius, + pixel_scales=0.1, + centre=centre, + ) + + adapt_data = aa.Array2D.ones( + shape_native=mask.shape_native, + pixel_scales=0.1, + ) + + kmeans = aa.image_mesh.Hilbert(pixels=16) + image_mesh = kmeans.image_plane_mesh_grid_from(mask=mask, adapt_data=adapt_data) + + points = np.asarray(image_mesh) + distances = np.sqrt( + (points[:, 0] - centre[0]) ** 2 + (points[:, 1] - centre[1]) ** 2 + ) + + assert np.all(distances <= radius + 1e-6) diff --git a/test_autoarray/mask/test_mask_2d.py b/test_autoarray/mask/test_mask_2d.py index 6bd9f4d51..5b75a60ac 100644 --- a/test_autoarray/mask/test_mask_2d.py +++ b/test_autoarray/mask/test_mask_2d.py @@ -705,3 +705,63 @@ def test__circular_radius__circular_mask__returns_radius_used_to_create_mask( ) assert mask.circular_radius == pytest.approx(radius, 1e-4) + + +@pytest.mark.parametrize( + "centre", + [ + (0.25, 0.0), + (0.0, 0.25), + (0.27, 0.13), + (-0.34, 0.41), + (0.5, 0.5), + ], +) +def test__is_circular__offset_centre_circular_mask__returns_true(centre): + mask = aa.Mask2D.circular( + shape_native=(20, 20), + radius=0.4, + pixel_scales=(0.1, 0.1), + centre=centre, + ) + + assert mask.is_circular == True + + +def test__is_circular__annular_mask__returns_false(): + mask = aa.Mask2D.circular_annular( + shape_native=(20, 20), + inner_radius=0.1, + outer_radius=0.4, + pixel_scales=(0.1, 0.1), + ) + + assert mask.is_circular == False + + +def test__is_circular__square_unmasked_region__returns_false(): + mask_array = np.full((10, 10), True) + mask_array[3:8, 3:8] = False + + mask = aa.Mask2D(mask=mask_array, pixel_scales=(1.0, 1.0)) + + assert mask.is_circular == False + + +@pytest.mark.parametrize( + "centre", + [ + (0.25, 0.0), + (0.0, 0.25), + (-0.34, 0.41), + ], +) +def test__circular_radius__offset_centre_circular_mask__returns_radius(centre): + mask = aa.Mask2D.circular( + shape_native=(40, 40), + radius=0.8, + pixel_scales=(0.1, 0.1), + centre=centre, + ) + + assert mask.circular_radius == pytest.approx(0.8, abs=0.1)