diff --git a/autoarray/__init__.py b/autoarray/__init__.py
index 326e4d10d..e1fa48a2e 100644
--- a/autoarray/__init__.py
+++ b/autoarray/__init__.py
@@ -64,7 +64,6 @@
 from .operators.contour import Grid2DContour
 from .layout.layout import Layout1D
 from .layout.layout import Layout2D
-from .preloads import Preloads
 from .structures.arrays.uniform_1d import Array1D
 from .structures.arrays.uniform_2d import Array2D
 from .structures.arrays.rgb import Array2DRGB
diff --git a/autoarray/config/general.yaml b/autoarray/config/general.yaml
index 257cae846..e0b0df0de 100644
--- a/autoarray/config/general.yaml
+++ b/autoarray/config/general.yaml
@@ -5,6 +5,7 @@ psf:
 inversion:
   check_reconstruction: true          # If True, the inversion's reconstruction is checked to ensure the solution of a meshs's mapper is not an invalid solution where the values are all the same.
   use_positive_only_solver: true      # If True, inversion's use a positive-only linear algebra solver by default, which is slower but prevents unphysical negative values in the reconstructed solutuion.
+  use_edge_zeroed_pixels : true       # If True, the edge pixels of a pixelization are set to zero, which prevents unphysical values in the reconstructed solution at the edge of the pixelization.
   no_regularization_add_to_curvature_diag_value : 1.0e-3 # The default value added to the curvature matrix's diagonal when regularization is not applied to a linear object, which prevents inversion's failing due to the matrix being singular.
   use_border_relocator: false          # If True, by default a pixelization's border is used to relocate all pixels outside its border to the border.
   reconstruction_vmax_factor: 0.5     # Plots of an Inversion's reconstruction use the reconstructed data's bright value multiplied by this factor.
diff --git a/autoarray/fixtures.py b/autoarray/fixtures.py
index fcba86a48..47598bbec 100644
--- a/autoarray/fixtures.py
+++ b/autoarray/fixtures.py
@@ -410,7 +410,7 @@ def make_delaunay_mapper_9_3x3():
         pixel_scales=1.0,
     )
 
-    mesh = aa.mesh.Delaunay()
+    mesh = aa.mesh.Delaunay(pixels=9)
 
     interpolator = mesh.interpolator_from(
         source_plane_data_grid=make_grid_2d_sub_2_7x7(),
@@ -443,7 +443,7 @@ def make_knn_mapper_9_3x3():
         pixel_scales=1.0,
     )
 
-    mesh = aa.mesh.KNearestNeighbor(split_neighbor_division=1)
+    mesh = aa.mesh.KNearestNeighbor(pixels=9, split_neighbor_division=1)
 
     interpolator = mesh.interpolator_from(
         source_plane_data_grid=make_grid_2d_sub_2_7x7(),
diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py
index af5cd189e..b45b1984a 100644
--- a/autoarray/inversion/inversion/abstract.py
+++ b/autoarray/inversion/inversion/abstract.py
@@ -12,7 +12,6 @@
 from autoarray.inversion.mappers.abstract import Mapper
 from autoarray.inversion.regularization.abstract import AbstractRegularization
 from autoarray.settings import Settings
-from autoarray.preloads import Preloads
 from autoarray.structures.arrays.uniform_2d import Array2D
 from autoarray.structures.grids.irregular_2d import Grid2DIrregular
 from autoarray.structures.visibilities import Visibilities
@@ -27,7 +26,6 @@ def __init__(
         dataset: Union[Imaging, Interferometer, DatasetInterface],
         linear_obj_list: List[LinearObj],
         settings: Settings = None,
-        preloads: Preloads = None,
         xp=np,
     ):
         """
@@ -74,8 +72,6 @@ def __init__(
 
         self.settings = settings or Settings()
 
-        self.preloads = preloads or Preloads()
-
         self.use_jax = xp is not np
 
     @property
@@ -234,9 +230,6 @@ def no_regularization_index_list(self) -> List[int]:
     @property
     def mapper_indices(self) -> np.ndarray:
 
-        if self.preloads.mapper_indices is not None:
-            return self.preloads.mapper_indices
-
         mapper_indices = []
 
         param_range_list = self.param_range_list_from(cls=Mapper)
@@ -386,6 +379,107 @@ def curvature_reg_matrix_reduced(self) -> Optional[np.ndarray]:
         # Zero rows and columns in the matrix we want to ignore
         return self.curvature_reg_matrix[ids_to_keep][:, ids_to_keep]
 
+    @cached_property
+    def zeroed_ids_to_keep(self):
+        """
+        Return the **positive global indices** of linear parameters that should be
+        kept (solved for) in the inversion, accounting for **zeroed pixel indices**
+        from one or more mappers.
+
+        ---------------------------------------------------------------------------
+        Parameter vector layout
+        ---------------------------------------------------------------------------
+        This method assumes the full linear parameter vector is ordered as:
+
+            [ non-pixel linear objects ][ mapper_0 pixels ][ mapper_1 pixels ] ... [ mapper_M pixels ]
+
+        where:
+
+        - *Non-pixel linear objects* include quantities such as analytic light
+          profiles, regularization amplitudes, etc.
+        - Each mapper contributes a contiguous block of pixel-based linear parameters.
+        - The concatenated pixel blocks occupy the **final** entries of the parameter
+          vector, with total length:
+
+              total_pixels = sum(mapper.mesh.pixels for mapper in mappers)
+
+        ---------------------------------------------------------------------------
+        Zeroed pixel convention
+        ---------------------------------------------------------------------------
+        For each mapper:
+
+        - `mapper.mesh.zeroed_pixels` must be a 1D array of **positive, mesh-local**
+          pixel indices in the range `[0, mapper.mesh.pixels - 1]`.
+        - These indices identify pixels that should be **excluded** from the linear
+          solve (e.g. edge pixels, masked regions, or padding pixels).
+        - Indexing is defined purely within the mapper’s own pixelization (e.g.
+          row-major flattening for rectangular meshes).
+
+        This method converts all mesh-local zeroed pixel indices into **global
+        parameter indices**, correctly offsetting for:
+          - the presence of non-pixel linear objects at the start of the vector
+          - the cumulative pixel counts of preceding mappers
+
+        ---------------------------------------------------------------------------
+        Backend and implementation details
+        ---------------------------------------------------------------------------
+        - The implementation is backend-agnostic and supports both NumPy and JAX via
+          `self._xp`.
+        - The returned indices are **positive global indices**, suitable for advanced
+          indexing of:
+              - `self.data_vector`
+              - `self.curvature_reg_matrix`
+        - When using JAX, this method avoids backend-incompatible operations and
+          preserves JIT compatibility under the same constraints as the rest of the
+          inversion pipeline.
+
+        Returns
+        -------
+        array-like
+            A 1D array of **positive global indices**, sorted in ascending order,
+            corresponding to linear parameters that should be kept in the inversion.
+        """
+
+        mapper_list = self.cls_list_from(cls=Mapper)
+
+        n_total = int(self.total_params)
+
+        pixels_per_mapper = [int(m.mesh.pixels) for m in mapper_list]
+        total_pixels = int(sum(pixels_per_mapper))
+
+        # Global start index of concatenated pixel block
+        pixel_start = n_total - total_pixels
+
+        # Total number of zeroed pixels across all mappers (Python int => static)
+        total_zeroed = int(sum(len(m.mesh.zeroed_pixels) for m in mapper_list))
+        n_keep = int(n_total - total_zeroed)
+
+        # Build global indices-to-zero across all mappers
+        zeros_global_list = []
+        offset = 0
+        for m, n_pix in zip(mapper_list, pixels_per_mapper):
+            zeros_local = self._xp.asarray(m.mesh.zeroed_pixels, dtype=self._xp.int32)
+            zeros_global_list.append(pixel_start + offset + zeros_local)
+            offset += n_pix
+
+        zeros_global = (
+            self._xp.concatenate(zeros_global_list)
+            if len(zeros_global_list) > 0
+            else self._xp.asarray([], dtype=self._xp.int32)
+        )
+
+        keep = self._xp.ones((n_total,), dtype=bool)
+
+        if self._xp is np:
+            keep[zeros_global] = False
+            keep_ids = self._xp.nonzero(keep)[0]
+
+        else:
+            keep = keep.at[zeros_global].set(False)
+            keep_ids = self._xp.nonzero(keep, size=n_keep)[0]
+
+        return keep_ids
+
     @cached_property
     def reconstruction(self) -> np.ndarray:
         """
@@ -405,16 +499,13 @@ def reconstruction(self) -> np.ndarray:
 
         if self.settings.use_positive_only_solver:
 
-            if self.preloads.source_pixel_zeroed_indices is not None:
-
-                # ids of values which are not zeroed and therefore kept in soluiton, which is computed in preloads.
-                ids_to_keep = self.preloads.source_pixel_zeroed_indices_to_keep
+            if self.settings.use_edge_zeroed_pixels and self.has(cls=Mapper):
 
                 # Use advanced indexing to select rows/columns
-                data_vector = self.data_vector[ids_to_keep]
-                curvature_reg_matrix = self.curvature_reg_matrix[ids_to_keep][
-                    :, ids_to_keep
-                ]
+                data_vector = self.data_vector[self.zeroed_ids_to_keep]
+                curvature_reg_matrix = self.curvature_reg_matrix[
+                    self.zeroed_ids_to_keep
+                ][:, self.zeroed_ids_to_keep]
 
                 # Perform reconstruction via fnnls
                 reconstruction_partial = (
@@ -431,11 +522,11 @@ def reconstruction(self) -> np.ndarray:
 
                 # Scatter the partial solution back to the full shape
                 if self._xp.__name__.startswith("jax"):
-                    reconstruction = reconstruction.at[ids_to_keep].set(
+                    reconstruction = reconstruction.at[self.zeroed_ids_to_keep].set(
                         reconstruction_partial
                     )
                 else:
-                    reconstruction[ids_to_keep] = reconstruction_partial
+                    reconstruction[self.zeroed_ids_to_keep] = reconstruction_partial
 
                 return reconstruction
 
diff --git a/autoarray/inversion/inversion/factory.py b/autoarray/inversion/inversion/factory.py
index e9955127a..01b5cd4d4 100644
--- a/autoarray/inversion/inversion/factory.py
+++ b/autoarray/inversion/inversion/factory.py
@@ -24,7 +24,6 @@
     InversionImagingSparse,
 )
 from autoarray.settings import Settings
-from autoarray.preloads import Preloads
 from autoarray.structures.arrays.uniform_2d import Array2D
 
 
@@ -32,7 +31,6 @@ def inversion_from(
     dataset: Union[Imaging, Interferometer, DatasetInterface],
     linear_obj_list: List[LinearObj],
     settings: Settings = None,
-    preloads: Preloads = None,
     xp=np,
 ):
     """
@@ -68,7 +66,6 @@ def inversion_from(
             dataset=dataset,
             linear_obj_list=linear_obj_list,
             settings=settings,
-            preloads=preloads,
             xp=xp,
         )
 
@@ -81,7 +78,6 @@ def inversion_imaging_from(
     dataset,
     linear_obj_list: List[LinearObj],
     settings: Settings = None,
-    preloads: Preloads = None,
     xp=np,
 ):
     """
@@ -133,7 +129,6 @@ def inversion_imaging_from(
                 dataset=dataset,
                 linear_obj_list=linear_obj_list,
                 settings=settings,
-                preloads=preloads,
                 xp=xp,
             )
 
@@ -141,7 +136,6 @@ def inversion_imaging_from(
             dataset=dataset,
             linear_obj_list=linear_obj_list,
             settings=settings,
-            preloads=preloads,
             xp=xp,
         )
 
@@ -149,7 +143,6 @@ def inversion_imaging_from(
         dataset=dataset,
         linear_obj_list=linear_obj_list,
         settings=settings,
-        preloads=preloads,
         xp=xp,
     )
 
diff --git a/autoarray/inversion/inversion/imaging/abstract.py b/autoarray/inversion/inversion/imaging/abstract.py
index 649382059..ae7dc8f44 100644
--- a/autoarray/inversion/inversion/imaging/abstract.py
+++ b/autoarray/inversion/inversion/imaging/abstract.py
@@ -8,7 +8,6 @@
 from autoarray.inversion.inversion.abstract import AbstractInversion
 from autoarray.inversion.linear_obj.linear_obj import LinearObj
 from autoarray.settings import Settings
-from autoarray.preloads import Preloads
 
 from autoarray.inversion.inversion.imaging import inversion_imaging_util
 
@@ -19,7 +18,6 @@ def __init__(
         dataset: Union[Imaging, DatasetInterface],
         linear_obj_list: List[LinearObj],
         settings: Settings = None,
-        preloads: Preloads = None,
         xp=np,
     ):
         """
@@ -67,7 +65,6 @@ def __init__(
             dataset=dataset,
             linear_obj_list=linear_obj_list,
             settings=settings,
-            preloads=preloads,
             xp=xp,
         )
 
diff --git a/autoarray/inversion/inversion/imaging/mapping.py b/autoarray/inversion/inversion/imaging/mapping.py
index 39ac0501e..b71a07b43 100644
--- a/autoarray/inversion/inversion/imaging/mapping.py
+++ b/autoarray/inversion/inversion/imaging/mapping.py
@@ -9,7 +9,6 @@
 from autoarray.inversion.linear_obj.linear_obj import LinearObj
 from autoarray.inversion.mappers.abstract import Mapper
 from autoarray.settings import Settings
-from autoarray.preloads import Preloads
 from autoarray.structures.arrays.uniform_2d import Array2D
 
 from autoarray.inversion.inversion import inversion_util
@@ -22,7 +21,6 @@ def __init__(
         dataset: Union[Imaging, DatasetInterface],
         linear_obj_list: List[LinearObj],
         settings: Settings = None,
-        preloads: Preloads = None,
         xp=np,
     ):
         """
@@ -49,7 +47,6 @@ def __init__(
             dataset=dataset,
             linear_obj_list=linear_obj_list,
             settings=settings,
-            preloads=preloads,
             xp=xp,
         )
 
diff --git a/autoarray/inversion/inversion/imaging/sparse.py b/autoarray/inversion/inversion/imaging/sparse.py
index ba32eb850..1025622df 100644
--- a/autoarray/inversion/inversion/imaging/sparse.py
+++ b/autoarray/inversion/inversion/imaging/sparse.py
@@ -10,7 +10,6 @@
 from autoarray.settings import Settings
 from autoarray.inversion.linear_obj.func_list import AbstractLinearObjFuncList
 from autoarray.inversion.mappers.abstract import Mapper
-from autoarray.preloads import Preloads
 from autoarray.structures.arrays.uniform_2d import Array2D
 
 from autoarray.inversion.inversion.imaging import inversion_imaging_util
@@ -22,7 +21,6 @@ def __init__(
         dataset: Union[Imaging, DatasetInterface],
         linear_obj_list: List[LinearObj],
         settings: Settings = None,
-        preloads: Preloads = None,
         xp=np,
     ):
         """
@@ -49,7 +47,6 @@ def __init__(
             dataset=dataset,
             linear_obj_list=linear_obj_list,
             settings=settings,
-            preloads=preloads,
             xp=xp,
         )
 
diff --git a/autoarray/inversion/inversion/imaging_numba/sparse.py b/autoarray/inversion/inversion/imaging_numba/sparse.py
index 31aabf5ff..4d8a4193d 100644
--- a/autoarray/inversion/inversion/imaging_numba/sparse.py
+++ b/autoarray/inversion/inversion/imaging_numba/sparse.py
@@ -10,7 +10,6 @@
 from autoarray.settings import Settings
 from autoarray.inversion.linear_obj.func_list import AbstractLinearObjFuncList
 from autoarray.inversion.mappers.abstract import Mapper
-from autoarray.preloads import Preloads
 from autoarray.structures.arrays.uniform_2d import Array2D
 
 from autoarray.inversion.inversion.imaging_numba import inversion_imaging_numba_util
@@ -22,7 +21,6 @@ def __init__(
         dataset: Union[Imaging, DatasetInterface],
         linear_obj_list: List[LinearObj],
         settings: Settings = None,
-        preloads: Preloads = None,
         xp=np,
     ):
         """
@@ -49,7 +47,6 @@ def __init__(
             dataset=dataset,
             linear_obj_list=linear_obj_list,
             settings=settings,
-            preloads=preloads,
             xp=xp,
         )
 
diff --git a/autoarray/inversion/mappers/abstract.py b/autoarray/inversion/mappers/abstract.py
index 706cfea95..d6d832b89 100644
--- a/autoarray/inversion/mappers/abstract.py
+++ b/autoarray/inversion/mappers/abstract.py
@@ -23,7 +23,6 @@ def __init__(
         regularization: Optional[AbstractRegularization] = None,
         settings: Settings = None,
         image_plane_mesh_grid=None,
-        preloads=None,
         xp=np,
     ):
         """
@@ -96,7 +95,7 @@ def __init__(
         self.interpolator = interpolator
 
         self.image_plane_mesh_grid = image_plane_mesh_grid
-        self.preloads = preloads
+
         self.settings = settings or Settings()
 
     @property
@@ -111,6 +110,10 @@ def pixels(self) -> int:
     def mask(self):
         return self.source_plane_data_grid.mask
 
+    @property
+    def mesh(self):
+        return self.interpolator.mesh
+
     @property
     def mesh_geometry(self):
         return self.interpolator.mesh_geometry
diff --git a/autoarray/inversion/mesh/interpolator/abstract.py b/autoarray/inversion/mesh/interpolator/abstract.py
index df8790365..7e59d8533 100644
--- a/autoarray/inversion/mesh/interpolator/abstract.py
+++ b/autoarray/inversion/mesh/interpolator/abstract.py
@@ -9,14 +9,12 @@ def __init__(
         mesh_grid,
         data_grid,
         adapt_data: np.ndarray = None,
-        preloads=None,
         xp=np,
     ):
         self.mesh = mesh
         self.mesh_grid = mesh_grid
         self.data_grid = data_grid
         self.adapt_data = adapt_data
-        self.preloads = preloads
         self.use_jax = xp is not np
 
     @property
diff --git a/autoarray/inversion/mesh/interpolator/delaunay.py b/autoarray/inversion/mesh/interpolator/delaunay.py
index deff4e981..63142d242 100644
--- a/autoarray/inversion/mesh/interpolator/delaunay.py
+++ b/autoarray/inversion/mesh/interpolator/delaunay.py
@@ -368,7 +368,6 @@ def __init__(
         mesh_grid,
         data_grid,
         adapt_data: np.ndarray = None,
-        preloads=None,
         xp=np,
     ):
         """
@@ -403,7 +402,6 @@ def __init__(
             mesh_grid=mesh_grid,
             data_grid=data_grid,
             adapt_data=adapt_data,
-            preloads=preloads,
             xp=xp,
         )
 
diff --git a/autoarray/inversion/mesh/interpolator/knn.py b/autoarray/inversion/mesh/interpolator/knn.py
index 584a96cc9..79d8d2556 100644
--- a/autoarray/inversion/mesh/interpolator/knn.py
+++ b/autoarray/inversion/mesh/interpolator/knn.py
@@ -217,7 +217,6 @@ def _mappings_sizes_weights_split(self):
             mesh=self.mesh,
             mesh_grid=self.mesh_grid,
             data_grid=split_points,
-            preloads=self.preloads,
             xp=self._xp,
         )
 
diff --git a/autoarray/inversion/mesh/interpolator/rectangular.py b/autoarray/inversion/mesh/interpolator/rectangular.py
index 2462bfbe6..e26c63148 100644
--- a/autoarray/inversion/mesh/interpolator/rectangular.py
+++ b/autoarray/inversion/mesh/interpolator/rectangular.py
@@ -236,7 +236,6 @@ def __init__(
         data_grid,
         mesh_weight_map,
         adapt_data: np.ndarray = None,
-        preloads=None,
         xp=np,
     ):
         """
@@ -271,7 +270,6 @@ def __init__(
             mesh_grid=mesh_grid,
             data_grid=data_grid,
             adapt_data=adapt_data,
-            preloads=preloads,
             xp=xp,
         )
         self.mesh_weight_map = mesh_weight_map
diff --git a/autoarray/inversion/mesh/interpolator/rectangular_uniform.py b/autoarray/inversion/mesh/interpolator/rectangular_uniform.py
index d3c7f7a3c..75d308d9b 100644
--- a/autoarray/inversion/mesh/interpolator/rectangular_uniform.py
+++ b/autoarray/inversion/mesh/interpolator/rectangular_uniform.py
@@ -108,7 +108,6 @@ def __init__(
         data_grid,
         adapt_data: np.ndarray = None,
         mesh_weight_map: np.ndarray = None,
-        preloads=None,
         xp=np,
     ):
         """
@@ -143,7 +142,6 @@ def __init__(
             mesh_grid=mesh_grid,
             data_grid=data_grid,
             adapt_data=adapt_data,
-            preloads=preloads,
             xp=xp,
         )
         self.mesh_weight_map = mesh_weight_map
diff --git a/autoarray/inversion/mesh/mesh/abstract.py b/autoarray/inversion/mesh/mesh/abstract.py
index 1b224cd9a..18120dc9f 100644
--- a/autoarray/inversion/mesh/mesh/abstract.py
+++ b/autoarray/inversion/mesh/mesh/abstract.py
@@ -99,7 +99,6 @@ def interpolator_from(
         source_plane_mesh_grid: Grid2DIrregular,
         border_relocator: Optional[BorderRelocator] = None,
         adapt_data: np.ndarray = None,
-        preloads=None,
         xp=np,
     ):
         raise NotImplementedError
diff --git a/autoarray/inversion/mesh/mesh/delaunay.py b/autoarray/inversion/mesh/mesh/delaunay.py
index 8da4bb671..0dc4f2ef7 100644
--- a/autoarray/inversion/mesh/mesh/delaunay.py
+++ b/autoarray/inversion/mesh/mesh/delaunay.py
@@ -8,32 +8,80 @@
 
 
 class Delaunay(AbstractMesh):
-    def __init__(self, areas_factor: float = 0.5):
+    def __init__(
+        self, pixels: int, zeroed_pixels: Optional[int] = 0, areas_factor: float = 0.5
+    ):
+        """
+        A Delaunay mesh composed of irregular triangular pixels used to reconstruct
+        a source on an unstructured grid.
+
+        The mesh consists of `pixels` vertices in the source plane, which are
+        connected via a Delaunay triangulation to form triangular elements.
+        Each vertex represents a linear parameter in the inversion.
+
+        Source-plane coordinates are interpolated onto this mesh using barycentric
+        interpolation within the enclosing triangle. For each coordinate, the three
+        vertices of the containing Delaunay triangle are identified and weighted
+        according to their barycentric distances, providing a smooth, piecewise-linear
+        reconstruction.
+
+        Zeroed pixels
+        -------------
+        The `zeroed_pixels` parameter specifies a number of mesh vertices that are
+        **excluded from the inversion**. These pixels are intended to correspond to
+        *edge or boundary vertices* of the Delaunay mesh.
+
+        Zeroing edge pixels helps to:
+          - stabilize the linear inversion,
+          - prevent poorly constrained boundary vertices from absorbing flux,
+          - reduce edge artefacts in the reconstructed source.
+
+        Zeroed pixels are always placed at the **end of the mesh parameter vector**
+        and are not solved for; their values are fixed to zero. Internally, the
+        inversion accounts for these excluded parameters when constructing and
+        solving the linear system.
+
+        Parameters
+        ----------
+        pixels : int
+            The number of active mesh vertices (linear parameters) used to represent
+            the source reconstruction.
+        areas_factor : float, optional
+            The barycentric area of Delaunay triangles is used to weight the regularization matrix.
+            This factor scales these areas, allowing for tuning of the regularization strength
+            based on triangle size.
+        zeroed_pixels : int, optional
+            The number of edge mesh vertices to exclude from the inversion. These
+            are appended to the end of the mesh and fixed to zero.
         """
-        An irregular mesh of Delaunay triangle pixels, which using linear barycentric interpolation are paired with
-        a 2D grid of (y,x) coordinates.
 
-        For a full description of how a mesh is paired with another grid,
-        see the :meth:`Pixelization API documentation <autoarray.inversion.pixelization.pixelization.Pixelization>`.
+        pixels = int(pixels) + zeroed_pixels
 
-        The Delaunay mesh represents pixels as an irregular 2D grid of Delaunay triangles.
+        super().__init__()
+        self.pixels = pixels
+        self.areas_factor = areas_factor
+        self._zeroed_pixels = zeroed_pixels
 
-        - ``image_plane_data_grid``: The observed data grid in the image-plane (which is paired with the mesh in
-          the source-plane).
-        - ``image_plane_mesh_grid``: The (y,x) mesh coordinates in the image-plane (which are the corners of Delaunay
-          triangles in the source-plane).
-        - ``source_plane_data_grid``: The observed data grid mapped to the source-plane after gravitational lensing.
-        - ``source_plane_mesh_grid``: The corner of each Delaunay triangle in the source-plane
-          (the ``image_plane_mesh_grid`` maps to this after gravitational lensing).
+    @property
+    def zeroed_pixels(self):
+        """
+        Return the **positive** mesh-local pixel indices to zero for a Delaunay mesh.
 
-        Each (y,x) coordinate in the ``source_plane_data_grid`` is paired with the three nearest Delaunay triangle
-        corners, using a weighted interpolation scheme.
+        For Delaunay meshes, `self.zeroed_pixels` is interpreted as a *count* of pixels
+        to be zeroed at the end of the pixel block. For example:
+            self.pixels = 780, self.zeroed_pixels = 30
+        returns indices 750..779.
 
-        Coordinates on the ``source_plane_data_grid`` are therefore given higher weights when paired with Delaunay
-        triangle corners they are a closer distance to.
+        Returns
+        -------
+        np.ndarray
+            1D array of positive pixel indices to zero.
         """
-        super().__init__()
-        self.areas_factor = areas_factor
+        if self._zeroed_pixels <= 0:
+            return np.array([], dtype=int)
+
+        start = self.pixels - self._zeroed_pixels
+        return np.arange(start, self.pixels, dtype=int)
 
     @property
     def skip_areas(self):
@@ -60,7 +108,6 @@ def interpolator_from(
         source_plane_mesh_grid: Grid2DIrregular,
         border_relocator: Optional[BorderRelocator] = None,
         adapt_data: np.ndarray = None,
-        preloads=None,
         xp=np,
     ):
         """
@@ -120,6 +167,5 @@ def interpolator_from(
             data_grid=relocated_grid,
             mesh_grid=relocated_mesh_grid,
             adapt_data=adapt_data,
-            preloads=preloads,
             xp=xp,
         )
diff --git a/autoarray/inversion/mesh/mesh/knn.py b/autoarray/inversion/mesh/mesh/knn.py
index 66dcdee24..929314a38 100644
--- a/autoarray/inversion/mesh/mesh/knn.py
+++ b/autoarray/inversion/mesh/mesh/knn.py
@@ -1,3 +1,5 @@
+from typing import Optional
+
 from autoarray.inversion.mesh.mesh.delaunay import Delaunay
 
 
@@ -5,18 +7,57 @@ class KNearestNeighbor(Delaunay):
 
     def __init__(
         self,
+        pixels: int,
+        zeroed_pixels: Optional[int] = 0,
         k_neighbors=10,
         radius_scale=1.5,
         areas_factor=0.5,
         split_neighbor_division=2,
     ):
+        """
+        A mesh that defines pixel connectivity using a k-nearest-neighbour
+        scheme rather than explicit triangle adjacency.
+
+        This mesh inherits the Delaunay geometry but does not use its interpolation scheme
+        but instead interpolates  by connecting each mesh vertex to its
+        `k_neighbors` nearest neighbouring vertices. These neighbour relationships are
+        used to impose smoothness constraints on the reconstructed source.
+
+        Neighbour connections may be further restricted using a distance-based criterion,
+        and optionally subdivided to improve stability for highly irregular meshes.
+
+        Parameters
+        ----------
+        pixels : int
+            The number of active mesh vertices (linear parameters) used to represent
+            the source reconstruction.
+        zeroed_pixels : int, optional
+            The number of edge mesh vertices to exclude from the inversion. These
+            boundary pixels are appended to the end of the parameter vector and
+            fixed to zero to reduce edge artefacts.
+        k_neighbors : int, optional
+            The number of nearest neighbours used to define connectivity for each
+            mesh vertex when constructing the regularization matrix.
+        radius_scale : float, optional
+            A multiplicative factor applied to the characteristic neighbour distance
+            that limits which neighbours are included. This prevents distant vertices
+            from contributing to regularization in sparsely sampled regions.
+        areas_factor : float, optional
+            The barycentric area of Delaunay triangles is used to weight the
+            regularization matrix. This factor scales these areas, allowing the
+            regularization strength to be tuned based on local triangle size.
+        split_neighbor_division : int, optional
+            Controls how neighbour connections are subdivided when forming the
+            regularization operator, improving numerical stability for irregular
+            point distributions.
+        """
 
         self.k_neighbors = k_neighbors
         self.radius_scale = radius_scale
         self.areas_factor = areas_factor
         self.split_neighbor_division = split_neighbor_division
 
-        super().__init__()
+        super().__init__(pixels=pixels, zeroed_pixels=zeroed_pixels)
 
     @property
     def skip_areas(self):
diff --git a/autoarray/inversion/mesh/mesh/rectangular_adapt_density.py b/autoarray/inversion/mesh/mesh/rectangular_adapt_density.py
index 61b0b1537..7d9006823 100644
--- a/autoarray/inversion/mesh/mesh/rectangular_adapt_density.py
+++ b/autoarray/inversion/mesh/mesh/rectangular_adapt_density.py
@@ -64,31 +64,53 @@ def overlay_grid_from(
 class RectangularAdaptDensity(AbstractMesh):
     def __init__(self, shape: Tuple[int, int] = (3, 3)):
         """
-        A uniform mesh of rectangular pixels, which without interpolation are paired with a 2D grid of (y,x)
-        coordinates.
+        A uniform rectangular mesh of pixels used to reconstruct a source on a
+        regular grid.
+
+        The mesh is defined by a 2D shape `(total_y_pixels, total_x_pixels)` and
+        is indexed in row-major order:
+
+            - Index 0 corresponds to the top-left pixel.
+            - Indices increase from left to right across each row,
+              and from top to bottom across rows.
+
+        Each source-plane coordinate is associated with the rectangular pixel
+        in which it lies. No interpolation is performed — every coordinate
+        contributes entirely to a single pixel.
+
+        Adaptive behaviour
+        ------------------
+        Although the rectangular mesh has a fixed, uniform geometry, it adapts
+        *implicitly* to the spatial density of the points it is paired with.
+        Regions of the source plane where many coordinates map onto the same
+        pixel receive stronger observational constraints, while sparsely
+        sampled regions are more weakly constrained.
+
+        In gravitational lensing applications, this naturally concentrates
+        information in regions of high magnification, where many image-plane
+        pixels map to a small area of the source plane. The inversion therefore
+        achieves higher effective resolution in these regions without requiring
+        explicit refinement of the mesh geometry.
+
+        Edge handling
+        -------------
+        Boundary (edge) pixels are automatically identified through the mesh
+        neighbour structure. These edge pixels may be internally excluded
+        (zeroed) during inversion to improve numerical stability and reduce
+        edge artefacts. This zeroing is determined by the mesh connectivity
+        and does not require manual specification of boundary indices.
 
-        For a full description of how a mesh is paired with another grid,
-        see the :meth:`Pixelization API documentation <autoarray.inversion.pixelization.pixelization.Pixelization>`.
-
-        The rectangular grid is uniform, has dimensions (total_y_pixels, total_x_pixels) and has indexing beginning
-        in the top-left corner and going rightwards and downwards.
-
-        A ``Pixelization`` using a ``RectangularAdaptDensity`` mesh has three grids associated with it:
-
-        - ``image_plane_data_grid``: The observed data grid in the image-plane (which is paired with the mesh in
-          the source-plane).
-        - ``source_plane_data_grid``: The observed data grid mapped to the source-plane after gravitational lensing.
-        - ``source_plane_mesh_grid``: The centres of each rectangular pixel.
-
-        It does not have a ``image_plane_mesh_grid`` because a rectangular pixelization is constructed by overlaying
-        a grid of rectangular over the `source_plane_data_grid`.
-
-        Each (y,x) coordinate in the `source_plane_data_grid` is associated with the rectangular pixelization pixel
-        it falls within. No interpolation is performed when making these associations.
         Parameters
         ----------
-        shape
-            The 2D dimensions of the rectangular grid of pixels (total_y_pixels, total_x_pixel).
+        shape : Tuple[int, int]
+            The 2D dimensions of the rectangular pixel grid
+            `(total_y_pixels, total_x_pixels)`.
+
+        Raises
+        ------
+        MeshException
+            If either dimension is less than 3, as a minimum of 3×3 pixels
+            is required to define interior and boundary structure.
         """
 
         if shape[0] <= 2 or shape[1] <= 2:
@@ -100,6 +122,32 @@ def __init__(self, shape: Tuple[int, int] = (3, 3)):
         self.pixels = self.shape[0] * self.shape[1]
         super().__init__()
 
+    @property
+    def zeroed_pixels(self):
+        """
+        Return the **positive** 1D pixel indices of the edge pixels in a rectangular mesh.
+
+        Indices are in row-major (C-order) flattened form for the rectangular pixel grid:
+            - 0 corresponds to the top-left pixel (row=0, col=0)
+            - indices increase across rows
+
+        These indices are defined purely within the rectangular mesh's pixel indexing
+        scheme (size = rows * cols) and are intended to be shifted / mapped to the full
+        inversion indexing inside the inversion logic.
+
+        Returns
+        -------
+        np.ndarray
+            A 1D array of positive indices corresponding to edge pixels.
+        """
+        from autoarray.inversion.mesh.mesh_geometry.rectangular import (
+            rectangular_edge_pixel_list_from,
+        )
+
+        edge_pixel_list = rectangular_edge_pixel_list_from(shape_native=self.shape)
+
+        return np.array(edge_pixel_list, dtype=int)
+
     @property
     def interpolator_cls(self):
         from autoarray.inversion.mesh.interpolator.rectangular import (
@@ -127,7 +175,6 @@ def interpolator_from(
         source_plane_mesh_grid: Grid2DIrregular,
         border_relocator: Optional[BorderRelocator] = None,
         adapt_data: np.ndarray = None,
-        preloads=None,
         xp=np,
     ):
         """
@@ -181,6 +228,5 @@ def interpolator_from(
             mesh_grid=Grid2DIrregular(mesh_grid),
             mesh_weight_map=mesh_weight_map,
             adapt_data=adapt_data,
-            preloads=preloads,
             xp=xp,
         )
diff --git a/autoarray/inversion/mesh/mesh/rectangular_adapt_image.py b/autoarray/inversion/mesh/mesh/rectangular_adapt_image.py
index 41a1badf1..0acd1a183 100644
--- a/autoarray/inversion/mesh/mesh/rectangular_adapt_image.py
+++ b/autoarray/inversion/mesh/mesh/rectangular_adapt_image.py
@@ -15,31 +15,63 @@ def __init__(
         weight_floor: float = 0.0,
     ):
         """
-        A uniform mesh of rectangular pixels, which without interpolation are paired with a 2D grid of (y,x)
-        coordinates.
+        A uniform rectangular mesh of pixels used to reconstruct a source on a
+        regular grid, with adaptive weighting driven by an external adapt image.
 
-        For a full description of how a mesh is paired with another grid,
-        see the :meth:`Pixelization API documentation <autoarray.inversion.pixelization.pixelization.Pixelization>`.
+        The mesh geometry is fixed and defined by a 2D shape
+        `(total_y_pixels, total_x_pixels)`. Pixels are indexed in row-major order:
 
-        The rectangular grid is uniform, has dimensions (total_y_pixels, total_x_pixels) and has indexing beginning
-        in the top-left corner and going rightwards and downwards.
+            - Index 0 corresponds to the top-left pixel.
+            - Indices increase left-to-right across rows and top-to-bottom
+              between rows.
 
-        A ``Pixelization`` using a ``RectangularAdaptDensity`` mesh has three grids associated with it:
+        Each source-plane coordinate is associated with the rectangular pixel
+        in which it lies. No interpolation is performed — every coordinate
+        contributes fully to a single pixel.
 
-        - ``image_plane_data_grid``: The observed data grid in the image-plane (which is paired with the mesh in
-          the source-plane).
-        - ``source_plane_data_grid``: The observed data grid mapped to the source-plane after gravitational lensing.
-        - ``source_plane_mesh_grid``: The centres of each rectangular pixel.
+        Adaptive behaviour (adapt image)
+        --------------------------------
+        Unlike a purely density-based rectangular mesh, this class adapts the
+        effective reconstruction using an *adapt image*. The adapt image provides
+        weights that emphasise specific regions of the source plane, typically
+        bright regions of a previously estimated reconstruction.
 
-        It does not have a ``image_plane_mesh_grid`` because a rectangular pixelization is constructed by overlaying
-        a grid of rectangular over the `source_plane_data_grid`.
+        Pixels corresponding to higher adapt-image intensity receive increased
+        weighting, allowing the inversion to prioritise reconstructing structure
+        in bright regions of the source. This leads to:
+
+          - improved resolution in high-signal regions,
+          - smoother behaviour in faint regions,
+          - reduced overfitting of noise in low-signal areas.
+
+        The weighting applied to each pixel is controlled by:
+
+          - `weight_power`: raises the adapt-image values to a power, increasing
+            or decreasing contrast between bright and faint regions.
+          - `weight_floor`: sets a minimum weight to prevent pixels in very faint
+            regions from becoming unconstrained.
+
+        This approach is particularly effective in strong gravitational lensing,
+        where the adapt image typically traces the intrinsic brightness
+        distribution of the source.
+
+        Edge handling
+        -------------
+        Boundary (edge) pixels are automatically identified via the mesh
+        neighbour structure and may be internally excluded (zeroed) during
+        inversion to improve numerical stability and reduce edge artefacts.
 
-        Each (y,x) coordinate in the `source_plane_data_grid` is associated with the rectangular pixelization pixel
-        it falls within. No interpolation is performed when making these associations.
         Parameters
         ----------
-        shape
-            The 2D dimensions of the rectangular grid of pixels (total_y_pixels, total_x_pixel).
+        shape : Tuple[int, int]
+            The 2D dimensions of the rectangular pixel grid
+            `(total_y_pixels, total_x_pixels)`.
+        weight_power : float, optional
+            Exponent applied to the adapt-image weights to control the strength
+            of adaptivity.
+        weight_floor : float, optional
+            Minimum weight applied to ensure numerical stability in low-intensity
+            regions.
         """
 
         super().__init__(shape=shape)
diff --git a/autoarray/inversion/mesh/mesh/rectangular_uniform.py b/autoarray/inversion/mesh/mesh/rectangular_uniform.py
index 6bf388987..d6c26d1e7 100644
--- a/autoarray/inversion/mesh/mesh/rectangular_uniform.py
+++ b/autoarray/inversion/mesh/mesh/rectangular_uniform.py
@@ -4,6 +4,39 @@
 
 
 class RectangularUniform(RectangularAdaptDensity):
+    """
+    A uniform rectangular mesh of pixels used to reconstruct a source on a
+    regular grid, with no adaptive weighting.
+
+    The mesh geometry is fixed and defined by a 2D shape
+    `(total_y_pixels, total_x_pixels)`. Pixels are indexed in row-major order:
+
+        - Index 0 corresponds to the top-left pixel.
+        - Indices increase left-to-right across rows and top-to-bottom
+          between rows.
+
+    Each source-plane coordinate is associated with the rectangular pixel
+    in which it lies. No interpolation is performed — every coordinate
+    contributes fully to a single pixel.
+
+    Uniform behaviour
+    -----------------
+    Unlike `RectangularAdaptDensity` and `RectangularAdaptImage`, this mesh
+    applies no adaptive weighting based on data density or an adapt image.
+    All pixels are treated equally in the reconstruction, and the effective
+    resolution is determined solely by the fixed mesh geometry and the
+    observational constraints.
+
+    This provides a simple and stable baseline reconstruction method,
+    particularly useful for controlled experiments or when adaptive
+    refinement is not required.
+
+    Edge handling
+    -------------
+    Boundary (edge) pixels are automatically identified via the mesh
+    neighbour structure and may be internally excluded (zeroed) during
+    inversion to improve numerical stability and reduce edge artefacts.
+    """
 
     @property
     def interpolator_cls(self):
diff --git a/autoarray/inversion/mesh/mesh_geometry/rectangular.py b/autoarray/inversion/mesh/mesh_geometry/rectangular.py
index 83bdda4db..c24d3f514 100644
--- a/autoarray/inversion/mesh/mesh_geometry/rectangular.py
+++ b/autoarray/inversion/mesh/mesh_geometry/rectangular.py
@@ -353,7 +353,7 @@ def rectangular_edges_from(shape_native, pixel_scales, xp=np):
 
 
 def rectangular_edge_pixel_list_from(
-    shape_native: Tuple[int, int], total_linear_light_profiles: int = 0
+    shape_native: Tuple[int, int],
 ) -> List[int]:
     """
     Returns a list of the 1D indices of all pixels on the edge of a rectangular pixelization,
@@ -382,12 +382,8 @@ def rectangular_edge_pixel_list_from(
     # Right column (excluding corners)
     right = (np.arange(1, rows - 1) + 1) * cols - 1
 
-    # Concatenate all edge indices
-    edge_pixel_indices = total_linear_light_profiles + np.concatenate(
-        [top, left, right, bottom]
-    )
+    edge_pixel_indices = np.concatenate([top, left, right, bottom])
 
-    # Sort and return
     return np.sort(edge_pixel_indices).tolist()
 
 
diff --git a/autoarray/preloads.py b/autoarray/preloads.py
deleted file mode 100644
index c72ad7323..000000000
--- a/autoarray/preloads.py
+++ /dev/null
@@ -1,120 +0,0 @@
-import logging
-
-import numpy as np
-
-logger = logging.getLogger(__name__)
-
-logger.setLevel(level="INFO")
-
-
-def mapper_indices_from(total_linear_light_profiles, total_mapper_pixels):
-
-    return np.arange(
-        total_linear_light_profiles,
-        total_linear_light_profiles + total_mapper_pixels,
-        dtype=int,
-    )
-
-
-class Preloads:
-
-    def __init__(
-        self,
-        mapper_indices: np.ndarray = None,
-        source_pixel_zeroed_indices: np.ndarray = None,
-        image_plane_mesh_grid_list: np.ndarray = None,
-        linear_light_profile_blurred_mapping_matrix=None,
-    ):
-        """
-        Stores preloaded arrays and matrices used during pixelized linear inversions, improving both performance
-        and compatibility with JAX.
-
-        Some arrays (e.g. `mapper_indices`) are required to be defined before sampling begins, because JAX demands
-        that input shapes remain static. These are used during each inversion to ensure consistent matrix shapes
-        for all likelihood evaluations.
-
-        Other arrays (e.g. parts of the curvature matrix) are preloaded purely to improve performance. In cases where
-        the source model is fixed (e.g. when fitting only the lens light), sections of the curvature matrix do not
-        change and can be reused, avoiding redundant computation.
-
-        Returns a list of image-plane mesh-grids, which are image-plane grids defining the centres of the pixels of
-        the pixelization's mesh (e.g. the centres rectangular mesh pixels).
-
-        The `image_mesh` attribute of the pixelization object defines whether the centre of each mesh pixel are
-        determined in the image-plane. When this is the case, the pixelization therefore has an image-plane mesh-grid,
-        which needs to be computed before the inversion is performed.
-
-        This function iterates over all galaxies with pixelizations, determines which pixelizations have an
-        `image_mesh` and for these pixelizations computes the image-plane mesh-grid.
-
-        It returns a list of all image-plane mesh-grids, which in the functions `interpolator_from` and `mapper_galaxy_dict`
-        are grouped into a `Mapper` object with other information required to perform the inversion using the
-        pixelization.
-
-        The order of this list is not important, because the `linear_obj_galaxy_dict` function associates each
-        mapper object (and therefore image-plane mesh-grid) with the galaxy it belongs to and is therefore used
-        elsewhere in the code (e.g. the fit module) to match inversion results to galaxies.
-
-        Certain image meshes adapt their pixels to the dataset, for example congregating the pixels to the brightest
-        regions of the image. This requires that `adapt_images` are used when setting up the image-plane mesh-grid.
-        This function uses the `adapt_images` attribute of the `GalaxiesToInversion` object pass these images and
-        raise an error if they are not present.
-
-        Parameters
-        ----------
-        mapper_indices
-            The integer indices of mapper pixels in the inversion. Used to extract reduced matrices (e.g.
-            `curvature_matrix_reduced`) that compute the pixelized inversion's log evidence term, where the indicies
-            are requirred to separate the rows and columns of matrices from linear light profiles.
-        source_pixel_zeroed_indices
-            Indices of source pixels that should be set to zero in the reconstruction. These typically correspond to
-            outer-edge source-plane regions with no image-plane mapping (e.g. outside a circular mask), helping
-            separate the lens light from the pixelized source model.
-        image_plane_mesh_grid
-            The (y,x) coordinates of the image-plane mesh grid used by pixelizations that start from pixels
-            being defined in the image-plane (e.g. overlaying a uniform grid of pixels on the image-plane, which
-            make up Delaunay triangles in the source-plane).
-        linear_light_profile_blurred_mapping_matrix
-            The evaluated images of the linear light profiles that make up the blurred mapping matrix component of the
-            inversion, with the other component being the pixelization's pixels. These are fixed when the lens light
-            is fixed to the maximum likelihood solution, allowing the blurred mapping matrix to be preloaded, but
-            the intensity values will still be solved for during the inversion.
-        """
-        self.mapper_indices = None
-        self.source_pixel_zeroed_indices = None
-        self.source_pixel_zeroed_indices_to_keep = None
-        self.linear_light_profile_blurred_mapping_matrix = None
-
-        if mapper_indices is not None:
-
-            self.mapper_indices = np.array(mapper_indices)
-
-        if source_pixel_zeroed_indices is not None:
-
-            self.source_pixel_zeroed_indices = np.array(source_pixel_zeroed_indices)
-
-            ids_zeros = np.array(source_pixel_zeroed_indices, dtype=int)
-
-            values_to_solve = np.ones(np.max(mapper_indices) + 1, dtype=bool)
-            values_to_solve[ids_zeros] = False
-
-            self.source_pixel_zeroed_indices_to_keep = np.where(values_to_solve)[0]
-
-        if image_plane_mesh_grid_list is not None:
-
-            self.image_plane_mesh_grid_list = []
-
-            for image_plane_mesh_grid in image_plane_mesh_grid_list:
-
-                if image_plane_mesh_grid is not None:
-                    self.image_plane_mesh_grid_list.append(
-                        np.array(image_plane_mesh_grid)
-                    )
-                else:
-                    self.image_plane_mesh_grid_list.append(None)
-
-        if linear_light_profile_blurred_mapping_matrix is not None:
-
-            self.linear_light_profile_blurred_mapping_matrix = np.array(
-                linear_light_profile_blurred_mapping_matrix
-            )
diff --git a/autoarray/settings.py b/autoarray/settings.py
index 8f33c99af..1533d77a6 100644
--- a/autoarray/settings.py
+++ b/autoarray/settings.py
@@ -12,6 +12,7 @@ def __init__(
         self,
         use_mixed_precision: bool = False,
         use_positive_only_solver: Optional[bool] = None,
+        use_edge_zeroed_pixels: Optional[bool] = None,
         use_border_relocator: Optional[bool] = None,
         no_regularization_add_to_curvature_diag_value: float = None,
     ):
@@ -41,6 +42,7 @@ def __init__(
         """
         self.use_mixed_precision = use_mixed_precision
         self._use_positive_only_solver = use_positive_only_solver
+        self._use_edge_zeroed_pixels = use_edge_zeroed_pixels
         self._use_border_relocator = use_border_relocator
         self._no_regularization_add_to_curvature_diag_value = (
             no_regularization_add_to_curvature_diag_value
@@ -53,6 +55,13 @@ def use_positive_only_solver(self):
 
         return self._use_positive_only_solver
 
+    @property
+    def use_edge_zeroed_pixels(self):
+        if self._use_edge_zeroed_pixels is None:
+            return conf.instance["general"]["inversion"]["use_edge_zeroed_pixels"]
+
+        return self._use_edge_zeroed_pixels
+
     @property
     def use_border_relocator(self):
         if self._use_border_relocator is None:
diff --git a/test_autoarray/config/general.yaml b/test_autoarray/config/general.yaml
index 5dfe968d2..c7b32a751 100644
--- a/test_autoarray/config/general.yaml
+++ b/test_autoarray/config/general.yaml
@@ -11,6 +11,7 @@ adapt:
 inversion:
   check_reconstruction: false       # If True, the inversion's reconstruction is checked to ensure the solution of a meshs's mapper is not an invalid solution where the values are all the same.
   use_positive_only_solver: false    # If True, inversion's use a positive-only linear algebra solver by default, which is slower but prevents unphysical negative values in the reconstructed solutuion.
+  use_edge_zeroed_pixels : false       # If True, the edge pixels of a pixelization are set to zero, which prevents unphysical values in the reconstructed solution at the edge of the pixelization.
   no_regularization_add_to_curvature_diag_value : 1.0e-8 # The default value added to the curvature matrix's diagonal when regularization is not applied to a linear object, which prevents inversion's failing due to the matrix being singular.
 numba:
   nopython: true
diff --git a/test_autoarray/inversion/inversion/test_abstract.py b/test_autoarray/inversion/inversion/test_abstract.py
index 51b3b21b2..e91f88810 100644
--- a/test_autoarray/inversion/inversion/test_abstract.py
+++ b/test_autoarray/inversion/inversion/test_abstract.py
@@ -159,8 +159,8 @@ def test__curvature_matrix_via_sparse_operator__includes_source_interpolation__i
 
     grid = aa.Grid2D.from_mask(mask=mask, over_sample_size=1)
 
-    mesh_0 = aa.mesh.Delaunay()
-    mesh_1 = aa.mesh.Delaunay()
+    mesh_0 = aa.mesh.Delaunay(pixels=9)
+    mesh_1 = aa.mesh.Delaunay(pixels=16)
 
     image_mesh_0 = aa.image_mesh.Overlay(shape=(3, 3))
     image_mesh_1 = aa.image_mesh.Overlay(shape=(4, 4))
@@ -559,7 +559,7 @@ def test__max_pixel_list_from_and_centre():
     mapper = aa.m.MockMapper(source_plane_mesh_grid=source_plane_mesh_grid)
 
     interpolator = aa.InterpolatorDelaunay(
-        mesh=aa.mesh.Delaunay(),
+        mesh=aa.mesh.Delaunay(pixels=4),
         mesh_grid=source_plane_mesh_grid,
         data_grid=None,
     )
@@ -597,7 +597,7 @@ def test__max_pixel_list_from__filter_neighbors():
     )
 
     mesh_geometry = aa.MeshGeometryDelaunay(
-        mesh=aa.mesh.Delaunay(),
+        mesh=aa.mesh.Delaunay(pixels=9),
         mesh_grid=source_plane_mesh_grid,
         data_grid=None,
     )
diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py
index a8b62b0f1..be751d07e 100644
--- a/test_autoarray/inversion/inversion/test_factory.py
+++ b/test_autoarray/inversion/inversion/test_factory.py
@@ -268,22 +268,22 @@ def test__inversion_imaging__via_regularizations(
     )
 
 
-def test__inversion_imaging__source_pixel_zeroed_indices(
+def test__inversion_imaging__zeroed_pixels(
     masked_imaging_7x7_no_blur,
     rectangular_mapper_7x7_3x3,
 ):
+
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[rectangular_mapper_7x7_3x3],
-        settings=aa.Settings(use_positive_only_solver=True),
-        preloads=aa.Preloads(
-            mapper_indices=range(0, 9), source_pixel_zeroed_indices=np.array([0])
+        settings=aa.Settings(
+            use_positive_only_solver=True, use_edge_zeroed_pixels=True
         ),
     )
 
     assert inversion.reconstruction.shape[0] == 9
     assert inversion.reconstruction[0] == 0.0
-    assert inversion.reconstruction[1] > 0.0
+    assert inversion.reconstruction[4] > 0.0
 
 
 def test__inversion_imaging__via_linear_obj_func_and_mapper(
diff --git a/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py b/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py
index aaa0e2980..6b775b006 100644
--- a/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py
+++ b/test_autoarray/inversion/pixelization/interpolator/test_delaunay.py
@@ -14,7 +14,7 @@ def test__scipy_delaunay__simplices(grid_2d_sub_1_7x7):
     )
 
     mesh_grid = aa.InterpolatorDelaunay(
-        mesh=aa.mesh.Delaunay(),
+        mesh=aa.mesh.Delaunay(pixels=6),
         mesh_grid=mesh_grid,
         data_grid=grid_2d_sub_1_7x7,
     )
@@ -34,7 +34,7 @@ def test__scipy_delaunay__split(grid_2d_sub_1_7x7):
     )
 
     mesh_grid = aa.InterpolatorDelaunay(
-        mesh=aa.mesh.Delaunay(),
+        mesh=aa.mesh.Delaunay(pixels=6),
         mesh_grid=mesh_grid,
         data_grid=grid_2d_sub_1_7x7,
     )
diff --git a/test_autoarray/inversion/pixelization/mappers/test_abstract.py b/test_autoarray/inversion/pixelization/mappers/test_abstract.py
index 0b0f258d3..58945cbba 100644
--- a/test_autoarray/inversion/pixelization/mappers/test_abstract.py
+++ b/test_autoarray/inversion/pixelization/mappers/test_abstract.py
@@ -139,7 +139,7 @@ def test__mapped_to_source_from(grid_2d_7x7):
         over_sample_size=1,
     )
 
-    mesh = aa.mesh.Delaunay()
+    mesh = aa.mesh.Delaunay(pixels=6)
 
     interpolator = mesh.interpolator_from(
         source_plane_data_grid=grid_2d_7x7,
diff --git a/test_autoarray/inversion/pixelization/mappers/test_delaunay.py b/test_autoarray/inversion/pixelization/mappers/test_delaunay.py
index f67586319..ebed69bfb 100644
--- a/test_autoarray/inversion/pixelization/mappers/test_delaunay.py
+++ b/test_autoarray/inversion/pixelization/mappers/test_delaunay.py
@@ -37,7 +37,7 @@ def test__pix_indexes_for_sub_slim_index__matches_util(grid_2d_sub_1_7x7):
         over_sample_size=1,
     )
 
-    mesh = aa.mesh.Delaunay()
+    mesh = aa.mesh.Delaunay(pixels=6)
 
     interpolator = mesh.interpolator_from(
         source_plane_data_grid=grid_2d_sub_1_7x7,
diff --git a/test_autoarray/inversion/pixelization/mappers/test_factory.py b/test_autoarray/inversion/pixelization/mappers/test_factory.py
index b767359c5..185f5b5bc 100644
--- a/test_autoarray/inversion/pixelization/mappers/test_factory.py
+++ b/test_autoarray/inversion/pixelization/mappers/test_factory.py
@@ -69,7 +69,7 @@ def test__delaunay_mapper():
     grid.over_sampled[0, 0] = -2.0
     grid.over_sampled[0, 1] = 2.0
 
-    mesh = aa.mesh.Delaunay()
+    mesh = aa.mesh.Delaunay(pixels=5)
     image_mesh = aa.image_mesh.Overlay(shape=(3, 3))
     image_plane_mesh_grid = image_mesh.image_plane_mesh_grid_from(
         mask=mask, adapt_data=None
diff --git a/test_autoarray/inversion/pixelization/mesh/test_abstract.py b/test_autoarray/inversion/pixelization/mesh/test_abstract.py
index f070a9b1f..cfbdbba2a 100644
--- a/test_autoarray/inversion/pixelization/mesh/test_abstract.py
+++ b/test_autoarray/inversion/pixelization/mesh/test_abstract.py
@@ -2,7 +2,7 @@
 
 
 def test__grid_is_relocated_via_border(grid_2d_7x7):
-    mesh = aa.mesh.Delaunay()
+    mesh = aa.mesh.Delaunay(pixels=9)
 
     mask = aa.Mask2D.circular(
         shape_native=(60, 60),
@@ -47,7 +47,7 @@ def test__grid_is_relocated_via_border(grid_2d_7x7):
     assert image_mesh[0, 0] != interpolator.mesh_grid[0, 0]
     assert interpolator.mesh_grid[0, 0] < 5.0
 
-    mesh = aa.mesh.Delaunay()
+    mesh = aa.mesh.Delaunay(pixels=9)
 
     border_relocator = aa.BorderRelocator(mask=mask, sub_size=1)
 
diff --git a/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py b/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py
index 3300bcb9d..0261421ef 100644
--- a/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py
+++ b/test_autoarray/inversion/pixelization/mesh_geometry/test_delaunay.py
@@ -14,7 +14,7 @@ def test__neighbors(grid_2d_sub_1_7x7):
     )
 
     mesh_geometry = aa.MeshGeometryDelaunay(
-        mesh=aa.mesh.Delaunay(),
+        mesh=aa.mesh.Delaunay(pixels=6),
         mesh_grid=mesh_grid,
         data_grid=grid_2d_sub_1_7x7,
     )
@@ -45,7 +45,7 @@ def test__voronoi_areas_via_delaunay_from(grid_2d_sub_1_7x7):
     )
 
     mesh = aa.MeshGeometryDelaunay(
-        mesh=aa.mesh.Delaunay(),
+        mesh=aa.mesh.Delaunay(pixels=6),
         mesh_grid=mesh_grid,
         data_grid=grid_2d_sub_1_7x7.over_sampled,
     )