Skip to content

Commit cef2871

Browse files
Jammy2211claude
authored andcommitted
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) <noreply@anthropic.com>
1 parent 3ddf368 commit cef2871

4 files changed

Lines changed: 160 additions & 19 deletions

File tree

autoarray/inversion/mesh/image_mesh/hilbert.py

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -118,10 +118,16 @@ def grid_hilbert_order_from(length, mask_radius):
118118
return x1d_hb, y1d_hb
119119

120120

121-
def image_and_grid_from(image, mask, mask_radius, pixel_scales, hilbert_length):
121+
def image_and_grid_from(
122+
image, mask, mask_radius, pixel_scales, hilbert_length, mask_centre=(0.0, 0.0)
123+
):
122124
"""
123125
This code will create a grid in Hilbert space-filling curve order and an interpolated hyper
124126
image associated to that grid.
127+
128+
The Hilbert curve is generated around the origin and filtered to a disc of radius ``mask_radius``, then translated
129+
by ``mask_centre`` so the sampling lands inside the mask's circle for offset-centre masks. For masks centred on
130+
the origin the translation is a no-op and behaviour is unchanged.
125131
"""
126132

127133
from scipy.interpolate import griddata
@@ -145,6 +151,7 @@ def image_and_grid_from(image, mask, mask_radius, pixel_scales, hilbert_length):
145151
grid_hb = np.stack((y1d_hb, x1d_hb), axis=-1)
146152
grid_hb_radius = np.sqrt(grid_hb[:, 0] ** 2.0 + grid_hb[:, 1] ** 2.0)
147153
new_grid = grid_hb[grid_hb_radius <= mask_radius]
154+
new_grid = new_grid + np.asarray(mask_centre)
148155

149156
new_img = griddata(
150157
points=grid,
@@ -271,6 +278,7 @@ def image_plane_mesh_grid_from(
271278
mask_radius=mask.circular_radius,
272279
pixel_scales=mask.pixel_scales,
273280
hilbert_length=193,
281+
mask_centre=mask.mask_centre,
274282
)
275283

276284
weight_map = self.weight_map_from(adapt_data=adapt_data_hb)

autoarray/mask/mask_2d.py

Lines changed: 63 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -938,10 +938,23 @@ def resized_from(self, new_shape, pad_value: int = 0.0) -> Mask2D:
938938
@property
939939
def is_circular(self) -> bool:
940940
"""
941-
Returns whether the mask is circular or not.
941+
Returns whether the unmasked region of the mask is a filled circular disc.
942942
943-
This is performed by taking the central row and column of the mask (based on the mask centre) and counting
944-
the number of unmasked pixels. If the number of unmasked pixels is the same, the mask is circular.
943+
The check is robust to circles whose centre is offset from the coordinate origin, including offsets that
944+
do not align with pixel boundaries. It rejects annular, square and other non-disc shapes.
945+
946+
The algorithm:
947+
948+
1) The bounding box of unmasked pixels must be square to within one pixel. Offset centres that fall between
949+
pixel centres can produce a one-pixel asymmetry in the bounding box, which is allowed; anything larger
950+
indicates a non-circular shape such as an ellipse. Tiny masks (bounding box <= 2 pixels) require an exactly
951+
square bounding box, since the one-pixel slack would otherwise admit 1x2 strips.
952+
2) The pixel at the geometric centre of the bounding box must be unmasked. This rejects annular masks whose
953+
inner hole is at least one pixel wide.
954+
3) The centre and radius of the unmasked region are inferred from the bounding box and pixel count, and a
955+
reference circular mask is built with those parameters. The input mask must match the reference within a
956+
small number of pixels (tolerance scales with mask area to absorb rim quantization). This rejects squares,
957+
crosses and tight annuli that slipped past the earlier checks.
945958
946959
This function does not support rectangular masks and an exception will be raised if the pixel scales in each
947960
direction are different.
@@ -955,22 +968,55 @@ def is_circular(self) -> bool:
955968
"""
956969
)
957970

958-
pixel_coordinates_2d = self.geometry.pixel_coordinates_2d_from(
959-
scaled_coordinates_2d=self.mask_centre
971+
where = np.where(np.invert(self.array))
972+
if where[0].size == 0:
973+
return False
974+
975+
y_min, y_max = int(where[0].min()), int(where[0].max())
976+
x_min, x_max = int(where[1].min()), int(where[1].max())
977+
y_extent = y_max - y_min + 1
978+
x_extent = x_max - x_min + 1
979+
980+
if abs(y_extent - x_extent) > 1:
981+
return False
982+
if max(y_extent, x_extent) <= 2 and y_extent != x_extent:
983+
return False
984+
985+
cy = (y_max + y_min) // 2
986+
cx = (x_max + x_min) // 2
987+
if bool(self.array[cy, cx]):
988+
return False
989+
990+
actual_area = int(where[0].size)
991+
inferred_radius_pix = np.sqrt(actual_area / np.pi)
992+
inferred_radius = inferred_radius_pix * self.pixel_scales[0]
993+
994+
y_centre_scaled = (
995+
0.5 * (self.shape_native[0] - 1) - 0.5 * (y_min + y_max)
996+
) * self.pixel_scales[0]
997+
x_centre_scaled = (
998+
0.5 * (x_min + x_max) - 0.5 * (self.shape_native[1] - 1)
999+
) * self.pixel_scales[1]
1000+
1001+
expected = mask_2d_util.mask_2d_circular_from(
1002+
shape_native=self.shape_native,
1003+
pixel_scales=self.pixel_scales,
1004+
radius=inferred_radius,
1005+
centre=(y_centre_scaled, x_centre_scaled),
9601006
)
9611007

962-
central_row_pixels = sum(np.invert(self[pixel_coordinates_2d[0], :]))
963-
central_column_pixels = sum(np.invert(self[:, pixel_coordinates_2d[1]]))
964-
965-
return central_row_pixels == central_column_pixels
1008+
diff_count = int(np.sum(self.array != expected))
1009+
tolerance = max(2, int(0.1 * actual_area))
1010+
return diff_count <= tolerance
9661011

9671012
@property
9681013
def circular_radius(self) -> float:
9691014
"""
9701015
Returns the radius in scaled units of a circular mask.
9711016
972-
This is performed by taking the central row of the mask (based on the mask centre) and counting the number of
973-
unmasked pixels. The radius is then half the number of unmasked pixels times the pixel scale.
1017+
The radius is computed from the bounding box of the unmasked region, taking the larger of the y and x extents
1018+
as the diameter in pixels. This is robust to offset centres that fall between pixel boundaries (where the
1019+
bounding box can be asymmetric by one pixel).
9741020
9751021
The mask is first checked that it is circular using the `is_circular` property, with an exception raised if
9761022
it is not.
@@ -987,15 +1033,14 @@ def circular_radius(self) -> float:
9871033
raise exc.MaskException(
9881034
"""
9891035
A circular radius can only be computed for a circular mask.
990-
1036+
9911037
The `is_circular` property of this mask has returned False, indicating the mask is not circular.
9921038
"""
9931039
)
9941040

995-
pixel_coordinates_2d = self.geometry.pixel_coordinates_2d_from(
996-
scaled_coordinates_2d=self.mask_centre
997-
)
998-
999-
central_row_pixels = sum(np.invert(self[pixel_coordinates_2d[0], :]))
1041+
where = np.where(np.invert(self.array))
1042+
y_extent = int(where[0].max() - where[0].min() + 1)
1043+
x_extent = int(where[1].max() - where[1].min() + 1)
1044+
diameter = max(y_extent, x_extent)
10001045

1001-
return central_row_pixels * self.pixel_scales[0] / 2.0
1046+
return diameter * self.pixel_scales[0] / 2.0

test_autoarray/inversion/pixelization/image_mesh/test_hilbert.py

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,4 @@
1+
import numpy as np
12
import pytest
23

34
import autoarray as aa
@@ -22,3 +23,30 @@ def test__image_plane_mesh_grid_from():
2223
[-1.02590674, -1.70984456],
2324
1.0e-4,
2425
)
26+
27+
28+
def test__image_plane_mesh_grid_from__offset_centre__points_inside_mask_circle():
29+
centre = (0.5, 0.3)
30+
radius = 0.5
31+
32+
mask = aa.Mask2D.circular(
33+
shape_native=(40, 40),
34+
radius=radius,
35+
pixel_scales=0.1,
36+
centre=centre,
37+
)
38+
39+
adapt_data = aa.Array2D.ones(
40+
shape_native=mask.shape_native,
41+
pixel_scales=0.1,
42+
)
43+
44+
kmeans = aa.image_mesh.Hilbert(pixels=16)
45+
image_mesh = kmeans.image_plane_mesh_grid_from(mask=mask, adapt_data=adapt_data)
46+
47+
points = np.asarray(image_mesh)
48+
distances = np.sqrt(
49+
(points[:, 0] - centre[0]) ** 2 + (points[:, 1] - centre[1]) ** 2
50+
)
51+
52+
assert np.all(distances <= radius + 1e-6)

test_autoarray/mask/test_mask_2d.py

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -705,3 +705,63 @@ def test__circular_radius__circular_mask__returns_radius_used_to_create_mask(
705705
)
706706

707707
assert mask.circular_radius == pytest.approx(radius, 1e-4)
708+
709+
710+
@pytest.mark.parametrize(
711+
"centre",
712+
[
713+
(0.25, 0.0),
714+
(0.0, 0.25),
715+
(0.27, 0.13),
716+
(-0.34, 0.41),
717+
(0.5, 0.5),
718+
],
719+
)
720+
def test__is_circular__offset_centre_circular_mask__returns_true(centre):
721+
mask = aa.Mask2D.circular(
722+
shape_native=(20, 20),
723+
radius=0.4,
724+
pixel_scales=(0.1, 0.1),
725+
centre=centre,
726+
)
727+
728+
assert mask.is_circular == True
729+
730+
731+
def test__is_circular__annular_mask__returns_false():
732+
mask = aa.Mask2D.circular_annular(
733+
shape_native=(20, 20),
734+
inner_radius=0.1,
735+
outer_radius=0.4,
736+
pixel_scales=(0.1, 0.1),
737+
)
738+
739+
assert mask.is_circular == False
740+
741+
742+
def test__is_circular__square_unmasked_region__returns_false():
743+
mask_array = np.full((10, 10), True)
744+
mask_array[3:8, 3:8] = False
745+
746+
mask = aa.Mask2D(mask=mask_array, pixel_scales=(1.0, 1.0))
747+
748+
assert mask.is_circular == False
749+
750+
751+
@pytest.mark.parametrize(
752+
"centre",
753+
[
754+
(0.25, 0.0),
755+
(0.0, 0.25),
756+
(-0.34, 0.41),
757+
],
758+
)
759+
def test__circular_radius__offset_centre_circular_mask__returns_radius(centre):
760+
mask = aa.Mask2D.circular(
761+
shape_native=(40, 40),
762+
radius=0.8,
763+
pixel_scales=(0.1, 0.1),
764+
centre=centre,
765+
)
766+
767+
assert mask.circular_radius == pytest.approx(0.8, abs=0.1)

0 commit comments

Comments
 (0)