Merge pull request #313 from PyAutoLabs/feature/hilbert-offset-centre-mask

fix(hilbert): support offset-centre circular masks

Author: Jammy2211

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)

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

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)

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)