diff --git a/autoarray/inversion/mappers/abstract.py b/autoarray/inversion/mappers/abstract.py
index b3c3e45fc..83a37d883 100644
--- a/autoarray/inversion/mappers/abstract.py
+++ b/autoarray/inversion/mappers/abstract.py
@@ -140,7 +140,11 @@ def source_plane_data_grid(self):
 
     @property
     def source_plane_mesh_grid(self):
-        return self.interpolator.mesh_grid
+        return self.interpolator.source_plane_mesh_grid
+
+    @property
+    def regularization_mesh_grid(self):
+        return self.interpolator.regularization_mesh_grid
 
     @property
     def adapt_data(self):
diff --git a/autoarray/inversion/mesh/interpolator/abstract.py b/autoarray/inversion/mesh/interpolator/abstract.py
index 2373a2acc..d193e4dde 100644
--- a/autoarray/inversion/mesh/interpolator/abstract.py
+++ b/autoarray/inversion/mesh/interpolator/abstract.py
@@ -30,6 +30,24 @@ def _xp(self):
             return jnp
         return np
 
+    @property
+    def source_plane_mesh_grid(self):
+        """
+        The source-plane (y,x) position of every mesh pixel, in the order of the reconstruction. For most meshes
+        this is the `mesh_grid` the interpolator was built from; the adaptive rectangular interpolator overrides it.
+        """
+        return self.mesh_grid
+
+    @property
+    def regularization_mesh_grid(self):
+        """
+        The (y,x) positions kernel regularizations measure pixel-to-pixel distances on. For most meshes these are
+        the source-plane pixel positions. For the adaptive rectangular meshes they are the regular lattice
+        `mesh_grid`, so distances are measured in the mesh's stretched (CDF) coordinates, see
+        `InterpolatorRectangular.source_plane_mesh_grid`.
+        """
+        return self.mesh_grid
+
     @property
     def _mappings_sizes_weights(self):
         raise NotImplementedError(
diff --git a/autoarray/inversion/mesh/interpolator/rectangular.py b/autoarray/inversion/mesh/interpolator/rectangular.py
index 527f9eefa..f87586162 100644
--- a/autoarray/inversion/mesh/interpolator/rectangular.py
+++ b/autoarray/inversion/mesh/interpolator/rectangular.py
@@ -478,6 +478,73 @@ def inv_full(U):
     return xp.abs(xp.outer(dy, dx).flatten())
 
 
+def adaptive_rectangular_node_centres_from(
+    source_grid_shape,
+    data_grid,
+    mesh_weight_map=None,
+    bandwidth: float = KERNEL_CDF_DEFAULT_BANDWIDTH,
+    n_knots: int = KERNEL_CDF_DEFAULT_KNOTS,
+    transform: str = "kernel",
+    xp=np,
+):
+    """
+    The source-plane (y,x) centre of every node of the adaptive rectangular mesh,
+    in the mapper's pixel order — the position where node ``(r, c)``'s
+    interpolation basis function peaks, so that reconstruction value ``k`` lives
+    at centre ``k``.
+
+    The mapper of ``adaptive_rectangular_mappings_weights_via_interpolation_from``
+    places the source plane in index space via ``g = (source_grid_size - 3) * U + 1``
+    with ``U`` the forward CDF transform, and its ``flatten`` puts y index-space
+    node ``k`` on row ``n - k`` and x node ``k`` on column ``k``. Node ``(r, c)``
+    is therefore at unit-square coordinate
+
+        U_y = (n - r - 1) / (n - 3),    U_x = (c - 1) / (n - 3),
+
+    mapped back through the SAME inverse CDF. The transform is separable, so each
+    axis is inverted once and the centres are their outer grid (row-major from
+    the top-left, like the uniform overlay this replaces).
+
+    Guard nodes (rows 0 and 1, columns 0 and ``n - 1``) have ``U`` outside
+    ``[0, 1]``; the inverse CDF clamps them to the edge of the data span, so they
+    sit on the border of the region the mapper samples (they carry zero
+    interpolation weight and zero area).
+
+    The mapper indexes a square mesh; for a non-square ``source_grid_shape`` each
+    axis uses its own pixel count, so there is still one centre per pixel.
+    """
+    n_y, n_x = source_grid_shape
+
+    mu = data_grid.mean(axis=0)
+    scale = data_grid.std(axis=0).min()
+    source_grid_scaled = (data_grid - mu) / scale
+
+    _, inv_transform = _transforms_from(
+        source_grid_scaled,
+        mesh_pixels=n_y,
+        mesh_weight_map=mesh_weight_map,
+        bandwidth=bandwidth,
+        n_knots=n_knots,
+        transform=transform,
+        xp=xp,
+    )
+
+    # the transform is per axis, so invert each axis's node coordinates in its own
+    # column (the other column is unused)
+    rows = xp.arange(n_y)
+    cols = xp.arange(n_x)
+    unit_y = (n_y - rows - 1) / (n_y - 3)
+    unit_x = (cols - 1) / (n_x - 3)
+
+    y_nodes = inv_transform(xp.stack([unit_y, xp.zeros(n_y)], axis=1))[:, 0]
+    x_nodes = inv_transform(xp.stack([xp.zeros(n_x), unit_x], axis=1))[:, 1]
+
+    y = xp.repeat(y_nodes * scale + mu[0], n_x)
+    x = xp.tile(x_nodes * scale + mu[1], n_y)
+
+    return xp.stack([y, x], axis=1)
+
+
 # Samples per axis of the 1D quadrature in `adaptive_rectangular_hat_areas_from`.
 # Against 40001 samples on a 12x12 mesh (PyAutoArray#659) the per-node error is
 # ~1e-6 on the smooth kernel CDF and ~6e-4 on the rank CDF's staircase over an
@@ -765,6 +832,23 @@ def mesh_geometry(self):
             xp=self._xp,
         )
 
+    @cached_property
+    def source_plane_mesh_grid(self):
+        """
+        The source-plane (y,x) centre of every mesh pixel, where the interpolant takes each reconstruction value
+        (``mesh_geometry.node_centres``). Use these to place, plot, export or interpolate the reconstruction.
+
+        ``mesh_grid`` is the uniform lattice before the adaptive CDF is applied. It is kept for the mesh
+        geometry's extent and is the ``regularization_mesh_grid``: kernel regularizations (``MaternKernel``,
+        ``GaussianKernel``, ...) and ``MaternSparse`` measure distances on it, i.e. in the mesh's stretched
+        coordinates, so a kernel's ``scale`` shrinks in the source plane where the mesh packs its pixels.
+        Distances on the true pixel centres gave a lower Bayesian evidence on a test lens, so the lattice is the
+        deliberate choice (PyAutoArray#667).
+        """
+        from autoarray.structures.grids.irregular_2d import Grid2DIrregular
+
+        return Grid2DIrregular(values=self.mesh_geometry.node_centres)
+
     @cached_property
     def _mappings_sizes_weights(self):
         mappings, weights = (
diff --git a/autoarray/inversion/mesh/mesh_geometry/rectangular.py b/autoarray/inversion/mesh/mesh_geometry/rectangular.py
index f323c0a33..9609258fc 100644
--- a/autoarray/inversion/mesh/mesh_geometry/rectangular.py
+++ b/autoarray/inversion/mesh/mesh_geometry/rectangular.py
@@ -496,6 +496,47 @@ def areas_for_magnification(self):
             xp=self._xp,
         )
 
+    @property
+    def node_centres(self):
+        """
+        The source-plane (y,x) centre of every pixel, in the order of the
+        reconstruction: the position at which the mapper's interpolant takes
+        each pixel's value.
+
+        For the adaptive meshes these are the nodes of the regular lattice in
+        CDF space mapped back through the mapper's inverse CDF (see
+        ``adaptive_rectangular_node_centres_from``), so they bunch up where the
+        mesh adapts. The uniform overlay ``mesh_grid`` is the lattice before
+        that mapping and is not where the values live. For
+        ``RectangularUniform`` the two are the same.
+        """
+        if self.transform == "uniform":
+            return self.mesh_grid.array
+
+        from autoarray.inversion.mesh.interpolator.rectangular import (
+            KERNEL_CDF_DEFAULT_BANDWIDTH,
+            KERNEL_CDF_DEFAULT_KNOTS,
+            adaptive_rectangular_node_centres_from,
+        )
+
+        return adaptive_rectangular_node_centres_from(
+            source_grid_shape=self.shape_native,
+            data_grid=self.data_grid.array,
+            mesh_weight_map=self.mesh_weight_map,
+            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
+            ),
+            transform=self.transform,
+            xp=self._xp,
+        )
+
     @property
     def edges_transformed(self):
         """
diff --git a/autoarray/inversion/mock/mock_mapper.py b/autoarray/inversion/mock/mock_mapper.py
index f73b0593b..e1eb7296e 100644
--- a/autoarray/inversion/mock/mock_mapper.py
+++ b/autoarray/inversion/mock/mock_mapper.py
@@ -51,6 +51,12 @@ def source_plane_mesh_grid(self):
             return super().source_plane_mesh_grid
         return self._source_plane_mesh_grid
 
+    @property
+    def regularization_mesh_grid(self):
+        if self._source_plane_mesh_grid is None:
+            return super().regularization_mesh_grid
+        return self._source_plane_mesh_grid
+
     @property
     def adapt_data(self):
         if self._adapt_data is None:
diff --git a/autoarray/inversion/regularization/abstract.py b/autoarray/inversion/regularization/abstract.py
index ae245e7da..8f30c7308 100644
--- a/autoarray/inversion/regularization/abstract.py
+++ b/autoarray/inversion/regularization/abstract.py
@@ -44,6 +44,22 @@ def validate_coefficient(coefficient, name: str = "coefficient"):
     )
 
 
+def regularization_mesh_grid_from(linear_obj):
+    """
+    The (y,x) positions a kernel regularization measures pixel-to-pixel distances on.
+
+    Mappers expose `regularization_mesh_grid` (the regular lattice for the adaptive rectangular meshes, see
+    `InterpolatorRectangular.source_plane_mesh_grid`). Other linear objects, including duck-typed ones in
+    downstream libraries, only define `source_plane_mesh_grid`, which is used for them.
+    """
+    grid = getattr(linear_obj, "regularization_mesh_grid", None)
+
+    if grid is None:
+        return linear_obj.source_plane_mesh_grid
+
+    return grid
+
+
 class AbstractRegularization:
     is_split_regularization = False
     """
diff --git a/autoarray/inversion/regularization/exponential_kernel.py b/autoarray/inversion/regularization/exponential_kernel.py
index 46203379e..ceb45e269 100644
--- a/autoarray/inversion/regularization/exponential_kernel.py
+++ b/autoarray/inversion/regularization/exponential_kernel.py
@@ -7,6 +7,7 @@
 
 from autoarray.inversion.regularization.abstract import AbstractRegularization
 from autoarray.inversion.regularization.abstract import validate_coefficient
+from autoarray.inversion.regularization.abstract import regularization_mesh_grid_from
 
 
 def exp_cov_matrix_from(
@@ -161,7 +162,7 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
 
         covariance_matrix = exp_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,
@@ -185,7 +186,7 @@ def log_det_regularization_matrix_term_from(
         """
         covariance_matrix = exp_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,
@@ -217,7 +218,7 @@ def regularization_term_from(
 
         covariance_matrix = exp_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,
diff --git a/autoarray/inversion/regularization/gaussian_kernel.py b/autoarray/inversion/regularization/gaussian_kernel.py
index 47beafef3..faa01663f 100644
--- a/autoarray/inversion/regularization/gaussian_kernel.py
+++ b/autoarray/inversion/regularization/gaussian_kernel.py
@@ -7,6 +7,7 @@
 
 from autoarray.inversion.regularization.abstract import AbstractRegularization
 from autoarray.inversion.regularization.abstract import validate_coefficient
+from autoarray.inversion.regularization.abstract import regularization_mesh_grid_from
 
 
 def gauss_cov_matrix_from(
@@ -150,7 +151,7 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
 
         covariance_matrix = gauss_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,
@@ -198,7 +199,7 @@ def log_det_regularization_matrix_term_from(
         """
         covariance_matrix = gauss_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,
@@ -236,7 +237,7 @@ def regularization_term_from(
 
         covariance_matrix = gauss_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
             xp=xp,
diff --git a/autoarray/inversion/regularization/matern_adapt_kernel.py b/autoarray/inversion/regularization/matern_adapt_kernel.py
index 8db0ce6d4..542ceac81 100644
--- a/autoarray/inversion/regularization/matern_adapt_kernel.py
+++ b/autoarray/inversion/regularization/matern_adapt_kernel.py
@@ -15,6 +15,7 @@
 )
 from autoarray.inversion.regularization.adapt import adapt_regularization_weights_from
 from autoarray.inversion.regularization.abstract import validate_coefficient
+from autoarray.inversion.regularization.abstract import regularization_mesh_grid_from
 
 
 class MaternAdaptKernel(MaternKernel):
@@ -146,7 +147,7 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
             linear_obj=linear_obj, xp=xp
         )
 
-        pixel_points = linear_obj.source_plane_mesh_grid.array
+        pixel_points = regularization_mesh_grid_from(linear_obj).array
 
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
@@ -178,7 +179,7 @@ def log_det_regularization_matrix_term_from(
 
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             nu=self.nu,
             weights=kernel_weights,
             jitter=self.jitter_value,
@@ -210,7 +211,7 @@ def regularization_term_from(
 
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             nu=self.nu,
             weights=kernel_weights,
             jitter=self.jitter_value,
diff --git a/autoarray/inversion/regularization/matern_kernel.py b/autoarray/inversion/regularization/matern_kernel.py
index 66b15f819..e70ff22c4 100644
--- a/autoarray/inversion/regularization/matern_kernel.py
+++ b/autoarray/inversion/regularization/matern_kernel.py
@@ -8,6 +8,7 @@
 
 from autoarray.inversion.regularization.abstract import AbstractRegularization
 from autoarray.inversion.regularization.abstract import validate_coefficient
+from autoarray.inversion.regularization.abstract import regularization_mesh_grid_from
 
 
 def kv_xp(v, z, xp=np):
@@ -367,7 +368,7 @@ def regularization_matrix_from(self, linear_obj: LinearObj, xp=np) -> np.ndarray
         """
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             nu=self.nu,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
@@ -390,7 +391,7 @@ def log_det_regularization_matrix_term_from(
         """
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             nu=self.nu,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
@@ -419,7 +420,7 @@ def regularization_term_from(
         """
         covariance_matrix = matern_cov_matrix_from(
             scale=self.scale,
-            pixel_points=linear_obj.source_plane_mesh_grid.array,
+            pixel_points=regularization_mesh_grid_from(linear_obj).array,
             nu=self.nu,
             jitter=self.jitter_value,
             jitter_relative=self.jitter_relative,
diff --git a/autoarray/inversion/regularization/matern_sparse.py b/autoarray/inversion/regularization/matern_sparse.py
index 2d7582483..ca260bdbb 100644
--- a/autoarray/inversion/regularization/matern_sparse.py
+++ b/autoarray/inversion/regularization/matern_sparse.py
@@ -24,6 +24,7 @@
 
 from autoarray.inversion.regularization.abstract import AbstractRegularization
 from autoarray.inversion.regularization.abstract import validate_coefficient
+from autoarray.inversion.regularization.abstract import regularization_mesh_grid_from
 
 
 def rectangular_laplacian_sparse_from(
@@ -33,8 +34,8 @@ def rectangular_laplacian_sparse_from(
     Returns the 5-point discrete Laplacian of a rectangular mesh with reflecting (zero-flux Neumann)
     edges, as a `scipy.sparse` CSR matrix of shape `(pixels, pixels)`.
 
-    Mesh pixels are indexed in row-major order, which is the order of a rectangular mesh's
-    `source_plane_mesh_grid` (built by `grid_2d_slim_via_shape_native_not_mask_from`): pixel
+    Mesh pixels are indexed in row-major order, which is the order of a rectangular mesh's lattice
+    (its `regularization_mesh_grid`, built by `grid_2d_slim_via_shape_native_not_mask_from`): pixel
     `k = i * total_x_pixels + j` is in row `i` (row 0 is the top row, the largest y) and column `j`
     (column 0 is the left column, the smallest x). Pixel `k` therefore neighbours `k +- 1` along x (within
     its row) and `k +- total_x_pixels` along y.
@@ -142,10 +143,11 @@ def __init__(
 
         **Meshes**: only rectangular meshes (`RectangularUniform`, `RectangularRTUAdaptDensity` and their
         subclasses) are supported, because the 5-point Laplacian needs the mesh's grid topology. The
-        Laplacian is built on the mapper's `source_plane_mesh_grid` (the uniform overlay grid, also the
-        grid the dense kernel measures distances on), whose pixels are in row-major order from the
-        top-left (see `rectangular_laplacian_sparse_from`). Other meshes (e.g. `Delaunay`) raise
-        `NotImplementedError`.
+        Laplacian is built on the mesh's `regularization_mesh_grid`, the regular lattice that the dense
+        kernels also measure distances on, whose pixels are in row-major order from the top-left (see
+        `rectangular_laplacian_sparse_from`). On an adaptive rectangular mesh that lattice is the mesh's
+        stretched coordinates, not the source-plane pixel centres (`source_plane_mesh_grid`). Other meshes
+        (e.g. `Delaunay`) raise `NotImplementedError`.
 
         **Smoothness**: only `nu = 1` is implemented. Higher integer `nu` is `c L^(nu + 1)`, a follow-up;
         non-integer `alpha` has no exact sparse precision.
@@ -243,7 +245,8 @@ def mesh_shape_and_pixel_scales_from(
     ) -> Tuple[Tuple[int, int], Tuple[float, float]]:
         """
         Returns the `(total_y_pixels, total_x_pixels)` shape and the `(y, x)` spacing of pixel centres of
-        a linear object's rectangular mesh, read from its `source_plane_mesh_grid`.
+        a linear object's rectangular mesh, read from its `regularization_mesh_grid` (the regular lattice, which
+        for the adaptive meshes is not the source-plane pixel centres `source_plane_mesh_grid`).
 
         Raises `NotImplementedError` if the linear object is not a mapper with a rectangular mesh.
         """
@@ -263,7 +266,7 @@ def mesh_shape_and_pixel_scales_from(
 
         shape_native = (int(mesh.shape[0]), int(mesh.shape[1]))
 
-        grid = np.asarray(linear_obj.source_plane_mesh_grid.array).reshape(
+        grid = np.asarray(regularization_mesh_grid_from(linear_obj).array).reshape(
             shape_native[0], shape_native[1], 2
         )
 
diff --git a/test_autoarray/inversion/pixelization/mesh_geometry/test_node_centres.py b/test_autoarray/inversion/pixelization/mesh_geometry/test_node_centres.py
new file mode 100644
index 000000000..f5ef82219
--- /dev/null
+++ b/test_autoarray/inversion/pixelization/mesh_geometry/test_node_centres.py
@@ -0,0 +1,217 @@
+"""
+The source-plane positions a rectangular mapper reports for its pixels
+(`mapper.source_plane_mesh_grid` = `mesh_geometry.node_centres`) must be where its
+interpolant takes each reconstruction value (PyAutoArray#667).
+
+For the adaptive meshes the interpolant lives on a regular lattice in CDF space, so
+the pixel centres are that lattice mapped back through the inverse CDF. The uniform
+overlay `mesh_grid` the interpolator is built from is the lattice before that mapping:
+using it as the pixel positions put the bright, densely packed central nodes on cells
+several times too large, so anything interpolating the reconstruction from those
+positions (the reconstruction CSV, `max_pixel_centre`, kernel regularization
+distances, the workspace's interpolated source) was wrong.
+"""
+
+import numpy as np
+import pytest
+from scipy.interpolate import griddata
+
+import autoarray as aa
+from autoarray.inversion.mesh.interpolator.rectangular import (
+    KERNEL_CDF_DEFAULT_BANDWIDTH,
+    KERNEL_CDF_DEFAULT_KNOTS,
+    adaptive_rectangular_mappings_weights_via_interpolation_from,
+)
+
+SHAPE = 24
+PIXEL_SCALE = 0.2
+SUB_SIZE = 2
+MESH_SHAPE = (12, 12)
+
+GAUSS_CENTRE = (0.2, -0.1)
+GAUSS_SIGMA = 0.5
+
+RECTANGULAR_ADAPTIVE = [
+    "RectangularBilinearAdaptDensity",
+    "RectangularBilinearAdaptImage",
+    "RectangularRTUAdaptDensity",
+    "RectangularRTUAdaptImage",
+]
+
+
+def _gaussian(grid):
+    return np.exp(
+        -((grid[:, 0] - GAUSS_CENTRE[0]) ** 2 + (grid[:, 1] - GAUSS_CENTRE[1]) ** 2)
+        / (2.0 * GAUSS_SIGMA**2)
+    )
+
+
+def _distort(points):
+    """
+    A smooth, lensing-like distortion of the uniform grid, so no two points share a
+    y or x coordinate. On an undistorted grid the coordinates are tied in whole rows
+    and columns, the rank CDF is a staircase and its inverse is ambiguous on the
+    steps (the documented worst case of the rank transform); a traced source-plane
+    grid has no such ties.
+    """
+    y, x = points[:, 0], points[:, 1]
+    return np.stack(
+        [y + 0.15 * np.sin(1.3 * x) + 0.05 * y**2, x + 0.15 * np.sin(1.1 * y)],
+        axis=1,
+    )
+
+
+def _mapper_from(name):
+    uniform = aa.Grid2D.uniform(
+        shape_native=(SHAPE, SHAPE),
+        pixel_scales=PIXEL_SCALE,
+        over_sample_size=SUB_SIZE,
+    )
+    grid = aa.Grid2D(
+        values=_distort(np.asarray(uniform.array)),
+        mask=uniform.mask,
+        over_sample_size=SUB_SIZE,
+        over_sampled=aa.Grid2DIrregular(
+            _distort(np.asarray(uniform.over_sampled.array))
+        ),
+    )
+    adapt_data = aa.Array2D(values=_gaussian(grid.array) + 0.01, mask=grid.mask)
+
+    interpolator = getattr(aa.mesh, name)(shape=MESH_SHAPE).interpolator_from(
+        source_plane_data_grid=grid,
+        source_plane_mesh_grid=None,
+        adapt_data=adapt_data,
+    )
+    return grid, aa.Mapper(interpolator=interpolator)
+
+
+def _interior():
+    n_y, n_x = MESH_SHAPE
+    interior = np.ones(MESH_SHAPE, dtype=bool)
+    interior[:2, :] = False
+    interior[:, 0] = False
+    interior[:, -1] = False
+    return interior.ravel()
+
+
+@pytest.mark.parametrize("name", RECTANGULAR_ADAPTIVE)
+def test__each_node_takes_its_full_weight_at_its_centre(name):
+    _, mapper = _mapper_from(name)
+
+    geometry = mapper.mesh_geometry
+    centres = np.asarray(mapper.source_plane_mesh_grid.array)
+
+    np.testing.assert_allclose(centres, np.asarray(geometry.node_centres))
+
+    mappings, weights = adaptive_rectangular_mappings_weights_via_interpolation_from(
+        source_grid_size=MESH_SHAPE[0],
+        data_grid=geometry.data_grid.array,
+        data_grid_over_sampled=centres,
+        mesh_weight_map=geometry.mesh_weight_map,
+        bandwidth=(
+            geometry.kernel_bandwidth
+            if geometry.kernel_bandwidth is not None
+            else KERNEL_CDF_DEFAULT_BANDWIDTH
+        ),
+        n_knots=(
+            geometry.kernel_knots
+            if geometry.kernel_knots is not None
+            else KERNEL_CDF_DEFAULT_KNOTS
+        ),
+        transform=geometry.transform,
+    )
+
+    pixels = np.arange(mapper.params)
+    self_weight = np.array([np.sum(weights[k][mappings[k] == k]) for k in pixels])
+
+    # The centres invert the mapper's forward CDF numerically, so a node's own weight
+    # at its centre is 1 up to that inversion: the kernel CDF's knot table (measured
+    # min 0.994 here) and the rank CDF's piecewise-linear interpolation between this
+    # coarse grid's 576 points (min 0.969).
+    tolerance = 0.04 if geometry.transform == "rank" else 0.01
+
+    assert self_weight[_interior()] == pytest.approx(1.0, abs=tolerance)
+
+
+@pytest.mark.parametrize("name", RECTANGULAR_ADAPTIVE)
+def test__centres_are_a_row_major_lattice_inside_the_data(name):
+    grid, mapper = _mapper_from(name)
+
+    centres = np.asarray(mapper.source_plane_mesh_grid.array).reshape(
+        MESH_SHAPE[0], MESH_SHAPE[1], 2
+    )
+
+    # y is constant along a row and decreases down the rows; x is constant down a
+    # column and increases along it
+    assert np.all(np.ptp(centres[:, :, 0], axis=1) == 0.0)
+    assert np.all(np.ptp(centres[:, :, 1], axis=0) == 0.0)
+    assert np.all(np.diff(centres[:, 0, 0]) <= 0.0)
+    assert np.all(np.diff(centres[0, :, 1]) >= 0.0)
+
+    # guard nodes are clamped to the data span, never outside it
+    data = np.asarray(grid.array)
+    assert np.all(centres[:, :, 0] >= data[:, 0].min() - 1.0e-12)
+    assert np.all(centres[:, :, 0] <= data[:, 0].max() + 1.0e-12)
+    assert np.all(centres[:, :, 1] >= data[:, 1].min() - 1.0e-12)
+    assert np.all(centres[:, :, 1] <= data[:, 1].max() + 1.0e-12)
+
+
+def _reconstruction_from(grid, mapper):
+    """A noise-free least-squares reconstruction of the Gaussian through the mapper."""
+    mapping_matrix = np.asarray(mapper.mapping_matrix_over_sampled)
+    data = _gaussian(np.asarray(grid.over_sampled.array))
+    curvature = mapping_matrix.T @ mapping_matrix
+    ridge = 1.0e-8 * np.mean(np.diag(curvature)) * np.eye(curvature.shape[0])
+    return np.linalg.solve(curvature + ridge, mapping_matrix.T @ data)
+
+
+def _interpolated_integral_from(points, reconstruction):
+    fine = aa.Grid2D.uniform(shape_native=(600, 600), pixel_scales=0.01)
+    interpolated = griddata(
+        points=points, values=reconstruction, xi=np.asarray(fine.array)
+    )
+    return np.nansum(interpolated) * 0.01**2
+
+
+@pytest.mark.parametrize("name", RECTANGULAR_ADAPTIVE)
+def test__interpolating_the_reconstruction_from_the_centres(name):
+    """
+    Linearly interpolating a reconstruction from the reported pixel positions must
+    put its flux in the right place. From the uniform overlay the integral was off by
+    a factor 0.59 to 3.7 on this setup.
+
+    On the density-adaptive meshes the integral then matches the analytic Gaussian to
+    ~2%. The image-adaptive meshes pack their nodes tightly onto the source, where
+    linear interpolation between centres in the source plane is a different
+    interpolant from the mapper's bilinear one in CDF space, so their integral is
+    only much closer than the overlay's (measured 1.30-1.64 against 2.9-3.7): the
+    mesh quadrature `areas_for_magnification` remains the way to integrate.
+    """
+    grid, mapper = _mapper_from(name)
+
+    reconstruction = _reconstruction_from(grid, mapper)
+    analytic = 2.0 * np.pi * GAUSS_SIGMA**2
+
+    centres = _interpolated_integral_from(
+        np.asarray(mapper.source_plane_mesh_grid.array), reconstruction
+    )
+    overlay = _interpolated_integral_from(
+        np.asarray(mapper.interpolator.mesh_grid.array), reconstruction
+    )
+
+    error_centres = abs(np.log(centres / analytic))
+    error_overlay = abs(np.log(overlay / analytic))
+
+    assert error_centres < 0.5 * error_overlay
+
+    if "AdaptDensity" in name:
+        assert centres / analytic == pytest.approx(1.0, abs=0.03)
+
+
+def test__rectangular_uniform__centres_are_the_overlay_grid():
+    _, mapper = _mapper_from("RectangularUniform")
+
+    np.testing.assert_array_equal(
+        np.asarray(mapper.source_plane_mesh_grid.array),
+        np.asarray(mapper.interpolator.mesh_grid.array),
+    )
diff --git a/test_autoarray/inversion/regularizations/test_matern_sparse.py b/test_autoarray/inversion/regularizations/test_matern_sparse.py
index 8dd8368d8..90ea4e711 100644
--- a/test_autoarray/inversion/regularizations/test_matern_sparse.py
+++ b/test_autoarray/inversion/regularizations/test_matern_sparse.py
@@ -41,7 +41,7 @@ def _laplacian_operator_by_hand_from(mapper, scale, nu=1.0):
     geometrically (same row and one x-spacing apart, or same column and one y-spacing apart) rather than
     from the row-major index arithmetic the class uses, so this also checks the mesh ordering.
     """
-    grid = np.asarray(mapper.source_plane_mesh_grid.array)
+    grid = np.asarray(mapper.regularization_mesh_grid.array)
     pixels = grid.shape[0]
 
     ys = np.unique(grid[:, 0])