Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 9 additions & 1 deletion autoarray/inversion/mesh/image_mesh/hilbert.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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,
Expand Down Expand Up @@ -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)
Expand Down
81 changes: 63 additions & 18 deletions autoarray/mask/mask_2d.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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.
Expand All @@ -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
28 changes: 28 additions & 0 deletions test_autoarray/inversion/pixelization/image_mesh/test_hilbert.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
import numpy as np
import pytest

import autoarray as aa
Expand All @@ -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)
60 changes: 60 additions & 0 deletions test_autoarray/mask/test_mask_2d.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Loading