From 20ed237a87e8cb80c30326a38c31b4a92fab497c Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Fri, 2 Oct 2026 11:32:51 +0100 Subject: [PATCH] fix: align transformed mesh areas with their cell edges --- .../mesh/mesh_geometry/rectangular.py | 36 +++--------- .../interpolator/test_rectangular.py | 34 ++++++++++- .../mesh_geometry/test_rectangular.py | 56 ++++++++++++++----- 3 files changed, 82 insertions(+), 44 deletions(-) diff --git a/autoarray/inversion/mesh/mesh_geometry/rectangular.py b/autoarray/inversion/mesh/mesh_geometry/rectangular.py index d11be49ad..1fb4775a6 100644 --- a/autoarray/inversion/mesh/mesh_geometry/rectangular.py +++ b/autoarray/inversion/mesh/mesh_geometry/rectangular.py @@ -437,36 +437,14 @@ def areas_transformed(self): The source-plane area of every pixel in the rectangular pixelization, transformed through the kernel-density CDF the mapper adapts with. - ``kernel_bandwidth`` / ``kernel_knots`` of None (e.g. the uniform - mesh's interpolator, which passes no kernel arguments) fall back to - the kernel defaults. + Cell boundaries are the mapper's node-midpoint edges, including + clamped guard cells. Reusing ``edges_transformed`` also keeps the + slim data grid aligned with its adapt-image weights. """ - from autoarray.inversion.mesh.interpolator.rectangular import ( - KERNEL_CDF_DEFAULT_BANDWIDTH, - KERNEL_CDF_DEFAULT_KNOTS, - adaptive_rectangular_areas_from, - ) - - bandwidth = ( - self.kernel_bandwidth - if self.kernel_bandwidth is not None - else KERNEL_CDF_DEFAULT_BANDWIDTH - ) - n_knots = ( - self.kernel_knots - if self.kernel_knots is not None - else KERNEL_CDF_DEFAULT_KNOTS - ) - - return adaptive_rectangular_areas_from( - source_grid_shape=self.shape_native, - data_grid=self.data_grid.over_sampled, - mesh_weight_map=self.mesh_weight_map, - bandwidth=bandwidth, - n_knots=n_knots, - transform=self.transform, - xp=self._xp, - ) + edges = self.edges_transformed + heights = self._xp.abs(self._xp.diff(edges[:, 0])) + widths = self._xp.abs(self._xp.diff(edges[:, 1])) + return self._xp.outer(heights, widths).ravel() @property def areas_for_magnification(self): diff --git a/test_autoarray/inversion/pixelization/interpolator/test_rectangular.py b/test_autoarray/inversion/pixelization/interpolator/test_rectangular.py index 485262aef..03a2d2d3e 100644 --- a/test_autoarray/inversion/pixelization/interpolator/test_rectangular.py +++ b/test_autoarray/inversion/pixelization/interpolator/test_rectangular.py @@ -471,7 +471,22 @@ def __getattr__(self, item): areas = geometry.areas_transformed assert areas.shape == (36,) assert np.all(np.isfinite(areas)) - assert np.all(areas > 0.0) + assert np.all(areas >= 0.0) + edges = geometry.edges_transformed + # For six nodes, y midpoints start at 11/6, 3/2, 7/6; + # x midpoints start at -1/2, -1/6 and end at 7/6, 3/2. + # Clamping to [0, 1] collapses two top rows and the side columns. + np.testing.assert_allclose(edges[:3, 0], np.max(data_grid.array[:, 0])) + np.testing.assert_allclose(edges[-1, 0], np.min(data_grid.array[:, 0])) + np.testing.assert_allclose(edges[:2, 1], np.min(data_grid.array[:, 1])) + np.testing.assert_allclose(edges[-2:, 1], np.max(data_grid.array[:, 1])) + areas_native = areas.reshape(6, 6) + np.testing.assert_array_equal(areas_native[:2], 0.0) + np.testing.assert_array_equal(areas_native[:, [0, -1]], 0.0) + assert np.all(areas_native[2:, 1:-1] > 0.0) + expected = np.outer(-np.diff(edges[:, 0]), np.diff(edges[:, 1])).ravel() + np.testing.assert_allclose(areas, expected) + assert areas.sum() == pytest.approx(np.prod(np.ptp(data_grid.array, axis=0))) def test__InterpolatorRectangular__rank_transform_via_property(): @@ -514,7 +529,22 @@ def __getattr__(self, item): areas = geometry.areas_transformed assert areas.shape == (36,) assert np.all(np.isfinite(areas)) - assert np.all(areas > 0.0) + assert np.all(areas >= 0.0) + edges = geometry.edges_transformed + # For six nodes, y midpoints start at 11/6, 3/2, 7/6; + # x midpoints start at -1/2, -1/6 and end at 7/6, 3/2. + # Clamping to [0, 1] collapses two top rows and the side columns. + np.testing.assert_allclose(edges[:3, 0], np.max(data_grid.array[:, 0])) + np.testing.assert_allclose(edges[-1, 0], np.min(data_grid.array[:, 0])) + np.testing.assert_allclose(edges[:2, 1], np.min(data_grid.array[:, 1])) + np.testing.assert_allclose(edges[-2:, 1], np.max(data_grid.array[:, 1])) + areas_native = areas.reshape(6, 6) + np.testing.assert_array_equal(areas_native[:2], 0.0) + np.testing.assert_array_equal(areas_native[:, [0, -1]], 0.0) + assert np.all(areas_native[2:, 1:-1] > 0.0) + expected = np.outer(-np.diff(edges[:, 0]), np.diff(edges[:, 1])).ravel() + np.testing.assert_allclose(areas, expected) + assert areas.sum() == pytest.approx(np.prod(np.ptp(data_grid.array, axis=0))) # --------------------------------------------------------------------------- diff --git a/test_autoarray/inversion/pixelization/mesh_geometry/test_rectangular.py b/test_autoarray/inversion/pixelization/mesh_geometry/test_rectangular.py index 370df69dd..075f6167a 100644 --- a/test_autoarray/inversion/pixelization/mesh_geometry/test_rectangular.py +++ b/test_autoarray/inversion/pixelization/mesh_geometry/test_rectangular.py @@ -159,23 +159,21 @@ def test__areas_transformed(mask_2d_7x7): ) areas = interpolator.mesh_geometry.areas_transformed + np.testing.assert_allclose( + interpolator.mesh_geometry.edges_transformed, + [[1.5, -1.5], [1.5, -1.5], [0.0, 0.0], [-1.5, 1.5]], + atol=1e-8, + ) assert np.all(np.isfinite(areas)) - assert np.all(areas > 0.0) - - # The unit square maps exactly onto the data bounding box (3.0 x 3.0), so - # the per-axis edge differences telescope and the areas sum to its area. + assert np.all(areas >= 0.0) + # Node-midpoint edges clamp the guard cells to the data span. The + # degenerate 3x3 mesh therefore covers its 3x3 bounding box with four cells. + np.testing.assert_allclose( + areas.reshape(3, 3), [[0.0, 0.0, 0.0], [0.0, 2.25, 2.25], [0.0, 2.25, 2.25]] + ) assert areas.sum() == pytest.approx(9.0, rel=1e-8) - # The input grid is 4-fold symmetric, so the transformed areas must be: - # equal corners, equal edge-midpoints, single centre value. - assert areas[0] == pytest.approx(areas[2], rel=1e-8) - assert areas[0] == pytest.approx(areas[6], rel=1e-8) - assert areas[0] == pytest.approx(areas[8], rel=1e-8) - assert areas[1] == pytest.approx(areas[3], rel=1e-8) - assert areas[1] == pytest.approx(areas[5], rel=1e-8) - assert areas[1] == pytest.approx(areas[7], rel=1e-8) - def test__edges_transformed(mask_2d_7x7): @@ -265,3 +263,35 @@ def test__edges_transformed__aligned_with_interpolation_node_convention(): # missed by ~0.4 in y. assert centroid_y == pytest.approx(test_point[0, 0], abs=0.1) assert centroid_x == pytest.approx(test_point[0, 1], abs=0.1) + + +@pytest.mark.parametrize("over_sample_size", [2, 1]) +def test__areas_transformed__adapt_image_matches_cell_edges(over_sample_size): + grid = aa.Grid2D.uniform( + shape_native=(5, 5), pixel_scales=0.4, over_sample_size=over_sample_size + ) + adapt_data = aa.Array2D.no_mask( + values=np.arange(1.0, 26.0).reshape(5, 5) ** 2, pixel_scales=0.4 + ) + mesh = aa.mesh.RectangularBilinearAdaptImage(shape=(6, 6)) + mapper = aa.Mapper( + interpolator=mesh.interpolator_from( + source_plane_data_grid=grid, + source_plane_mesh_grid=None, + adapt_data=adapt_data, + ) + ) + geometry = mapper.mesh_geometry + assert len(geometry.mesh_weight_map) == len(grid.array) + assert len(grid.over_sampled) == len(grid.array) * over_sample_size**2 + assert np.ptp(geometry.mesh_weight_map) > 0.0 + + edges = geometry.edges_transformed + areas = geometry.areas_transformed + # Rows descend in y; columns ascend in x. Preserve the mapper's row order. + expected = np.outer(-np.diff(edges[:, 0]), np.diff(edges[:, 1])).ravel() + assert areas.shape == (36,) + assert np.all(np.isfinite(areas)) + assert np.all(areas >= 0.0) + np.testing.assert_allclose(areas, expected, rtol=1e-12, atol=1e-12) + np.testing.assert_array_equal(geometry.areas_for_magnification, areas)