diff --git a/autoarray/config/general.yaml b/autoarray/config/general.yaml
index e4e42d3b7..6f17e06a3 100644
--- a/autoarray/config/general.yaml
+++ b/autoarray/config/general.yaml
@@ -1,6 +1,5 @@
 inversion:
   reconstruction_vmax_factor: 0.5   # Plots of an Inversion's reconstruction use the reconstructed data's bright value multiplied by this factor.
-  check_solution_default: true      # When an Inversion checks solutions in case they are numerically unstable by default or not (e.g. InversionSettings.check_solution).
 numba:
   use_numba: true
   cache: false
diff --git a/autoarray/config/visualize/general.yaml b/autoarray/config/visualize/general.yaml
index 4967b92b8..63ac89ea1 100644
--- a/autoarray/config/visualize/general.yaml
+++ b/autoarray/config/visualize/general.yaml
@@ -4,6 +4,5 @@ general:
   zoom_around_mask: true            # If True, plots of data structures with a mask automatically zoom in the masked region.
 inversion:
   reconstruction_vmax_factor: 0.5   # Plots of an Inversion's reconstruction use the reconstructed data's bright value multiplied by this factor.
-  check_solution_default: true      # When an Inversion checks solutions in case they are numerically unstable by default or not (e.g. InversionSettings.check_solution).
 units:
   in_kpc: false                     # If True, plots that are normally in arc-seconds are instead plotted using kiloparsecs.
\ No newline at end of file
diff --git a/autoarray/inversion/inversion/abstract.py b/autoarray/inversion/inversion/abstract.py
index 73e01fc94..6e589a540 100644
--- a/autoarray/inversion/inversion/abstract.py
+++ b/autoarray/inversion/inversion/abstract.py
@@ -1,3 +1,5 @@
+import copy
+
 import numpy as np
 from scipy.linalg import block_diag
 from scipy.sparse import csc_matrix
@@ -223,6 +225,33 @@ def all_linear_obj_have_regularization(self) -> bool:
             list(filter(None, self.regularization_list))
         )
 
+    @property
+    def mapper_edge_pixel_list(self) -> List[int]:
+        """
+        Returns the edge pixels of all mappers in the inversion.
+
+        This uses the `edge_pixel_list` property of the `Mesh` of the `Mapper` class, and updates their values to
+        correspond to the indexing of the overall inversion's `curvature_matrix`.
+
+        This is used to regulareze the edge pixels of the inversion's `reconstruction` or remove them from the
+        inversion procedure entirely (e.g. make these values of these edge pixels zero).
+
+        Returns
+        -------
+        A list of the edge pixels of all mappers in the inversion, where the values are updated to correspond to the
+        indexing of the overall inversion's `curvature_matrix`.
+        """
+        mapper_edge_pixel_list = []
+
+        param_range_list = self.param_range_list_from(cls=LinearObj)
+
+        for param_range, linear_obj in zip(param_range_list, self.linear_obj_list):
+            if isinstance(linear_obj, AbstractMapper):
+                for edge_pixel in linear_obj.edge_pixel_list:
+                    mapper_edge_pixel_list.append(edge_pixel + param_range[0])
+
+        return mapper_edge_pixel_list
+
     @property
     def total_regularizations(self) -> int:
         return sum(
@@ -323,6 +352,9 @@ def regularization_matrix(self) -> Optional[np.ndarray]:
         For multiple mappers, the regularization matrix is computed as the block diagonal of each individual mapper.
         The scipy function `block_diag` has an overhead associated with it and if there is only one mapper and
         regularization it is bypassed.
+
+        If the `settings.force_edge_pixels_to_zeros` is `True`, the edge pixels of each mapper in the inversion
+        are regularized so high their value is forced to zero.
         """
         if self.preloads.regularization_matrix is not None:
             return self.preloads.regularization_matrix
@@ -375,8 +407,6 @@ def curvature_reg_matrix(self) -> np.ndarray:
         if not self.has(cls=AbstractRegularization):
             return self.curvature_matrix
 
-        # TODO : Done for speed, instead check if there is one regularization
-
         if len(self.regularization_list) == 1:
 
             curvature_matrix = self.curvature_matrix
@@ -388,6 +418,19 @@ def curvature_reg_matrix(self) -> np.ndarray:
 
         return np.add(self.curvature_matrix, self.regularization_matrix)
 
+    @cached_property
+    def curvature_reg_matrix_solver(self):
+
+        if self.settings.force_edge_pixels_to_zeros:
+
+            curvature_reg_matrix_solver = copy.copy(self.curvature_reg_matrix)
+
+            curvature_reg_matrix_solver[:, self.mapper_edge_pixel_list] = 0.0
+
+            return curvature_reg_matrix_solver
+
+        return self.curvature_reg_matrix
+
     @cached_property
     @profile_func
     def curvature_reg_matrix_reduced(self) -> np.ndarray:
@@ -411,36 +454,6 @@ def curvature_reg_matrix_reduced(self) -> np.ndarray:
 
         return curvature_reg_matrix
 
-    @cached_property
-    @profile_func
-    def curvature_reg_matrix_cholesky(self) -> np.ndarray:
-        """
-        Performs a Cholesky decomposition of the `curvature_reg_matrix`, the result of which is used to solve the
-        linear system of equations of the `Inversion`.
-
-        The method `np.linalg.solve` is faster to do this, but the Cholesky decomposition is used later in the code
-        to speed up the calculation of `log_det_curvature_reg_matrix_term`.
-        """
-        try:
-            return np.linalg.cholesky(self.curvature_reg_matrix)
-        except np.linalg.LinAlgError:
-            raise exc.InversionException()
-
-    @cached_property
-    @profile_func
-    def curvature_reg_matrix_reduced_cholesky(self) -> np.ndarray:
-        """
-        Performs a Cholesky decomposition of the `curvature_reg_matrix`, the result of which is used to solve the
-        linear system of equations of the `Inversion`.
-
-        The method `np.linalg.solve` is faster to do this, but the Cholesky decomposition is used later in the code
-        to speed up the calculation of `log_det_curvature_reg_matrix_term`.
-        """
-        try:
-            return np.linalg.cholesky(self.curvature_reg_matrix_reduced)
-        except np.linalg.LinAlgError:
-            raise exc.InversionException()
-
     @cached_property
     @profile_func
     def reconstruction(self) -> np.ndarray:
@@ -453,12 +466,13 @@ def reconstruction(self) -> np.ndarray:
         if not self.settings.use_positive_only_solver:
             return inversion_util.reconstruction_positive_negative_from(
                 data_vector=self.data_vector,
-                curvature_reg_matrix_cholesky=self.curvature_reg_matrix_cholesky,
+                curvature_reg_matrix=self.curvature_reg_matrix_solver,
                 settings=self.settings,
             )
+
         return inversion_util.reconstruction_positive_only_from(
             data_vector=self.data_vector,
-            curvature_reg_matrix=self.curvature_reg_matrix,
+            curvature_reg_matrix=self.curvature_reg_matrix_solver,
             settings=self.settings,
         )
 
@@ -610,7 +624,12 @@ def log_det_curvature_reg_matrix_term(self) -> float:
         if not self.has(cls=AbstractRegularization):
             return 0.0
 
-        return 2.0 * np.sum(np.log(np.diag(self.curvature_reg_matrix_reduced_cholesky)))
+        try:
+            return 2.0 * np.sum(
+                np.log(np.diag(np.linalg.cholesky(self.curvature_reg_matrix_reduced)))
+            )
+        except np.linalg.LinAlgError:
+            raise exc.InversionException()
 
     @cached_property
     @profile_func
diff --git a/autoarray/inversion/inversion/factory.py b/autoarray/inversion/inversion/factory.py
index 64a506b2b..de4845b7e 100644
--- a/autoarray/inversion/inversion/factory.py
+++ b/autoarray/inversion/inversion/factory.py
@@ -227,7 +227,7 @@ def inversion_imaging_unpacked_from(
     An `Inversion` whose type is determined by the input `dataset` and `settings`.
     """
 
-    if all(
+    if any(
         isinstance(linear_obj, AbstractLinearObjFuncList)
         for linear_obj in linear_obj_list
     ):
diff --git a/autoarray/inversion/inversion/inversion_util.py b/autoarray/inversion/inversion/inversion_util.py
index 4d1a1f464..d277c219c 100644
--- a/autoarray/inversion/inversion/inversion_util.py
+++ b/autoarray/inversion/inversion/inversion_util.py
@@ -257,7 +257,7 @@ def mapped_reconstructed_data_via_mapping_matrix_from(
 
 def reconstruction_positive_negative_from(
     data_vector: np.ndarray,
-    curvature_reg_matrix_cholesky: np.ndarray,
+    curvature_reg_matrix: np.ndarray,
     settings: SettingsInversion = SettingsInversion(),
 ):
     """
@@ -278,8 +278,8 @@ def reconstruction_positive_negative_from(
     ----------
     data_vector
         The `data_vector` D which is solved for.
-    curvature_reg_matrix_cholesky
-        The cholesky decomposition of the sum of the curvature and regularization matrices.
+    curvature_reg_matrix
+        The sum of the curvature and regularization matrices.
     settings
         Controls the settings of the inversion, for this function where the solution is checked to not be all
         the same values.
@@ -290,14 +290,10 @@ def reconstruction_positive_negative_from(
         The curvature_matrix plus regularization matrix, overwriting the curvature_matrix in memory.
     """
     try:
-        reconstruction = cho_solve((curvature_reg_matrix_cholesky, True), data_vector)
+        reconstruction = np.linalg.solve(curvature_reg_matrix, data_vector)
     except np.linalg.LinAlgError as e:
         raise exc.InversionException() from e
 
-    reconstruction_check_solution(
-        data_vector=data_vector, reconstruction=reconstruction, settings=settings
-    )
-
     return reconstruction
 
 
@@ -324,7 +320,7 @@ def reconstruction_positive_only_from(
     ----------
     data_vector
         The `data_vector` D which is solved for.
-    curvature_reg_matrix_cholesky
+    curvature_reg_matrix
         The sum of the curvature and regularization matrices.
     settings
         Controls the settings of the inversion, for this function where the solution is checked to not be all
@@ -347,48 +343,9 @@ def reconstruction_positive_only_from(
     except (RuntimeError, np.linalg.LinAlgError) as e:
         raise exc.InversionException() from e
 
-    reconstruction_check_solution(
-        data_vector=data_vector, reconstruction=reconstruction, settings=settings
-    )
-
     return reconstruction
 
 
-def reconstruction_check_solution(
-    data_vector: np.ndarray,
-    reconstruction: np.ndarray,
-    settings: SettingsInversion = SettingsInversion(),
-):
-    """
-    Check that the solution after reconstructing an inversion is not one where all reconstructed values go to the same
-    value.
-
-    These solutions occur when the linear system is numerically unstable, and can lead to high `log_evidence` values
-    due to numerically issues not currently fully understood.
-
-    Parameters
-    ----------
-    data_vector
-        The `data_vector` D which is solved for.
-    reconstruction
-        The solution of the linear system whose solution is checked.
-    settings
-        Controls the settings of the inversion, for this function where the solution is checked to not be all
-        the same values.
-    """
-    if settings.check_solution and len(reconstruction) > 1:
-        if np.isclose(a=reconstruction[0], b=reconstruction[1], atol=1e-6).all():
-            raise exc.InversionException()
-
-    index = int(data_vector.shape[0] / 2)
-
-    if settings.check_solution and len(reconstruction) > 1:
-        if np.isclose(
-            a=reconstruction[index], b=reconstruction[index + 1], atol=1e-6
-        ).all():
-            raise exc.InversionException()
-
-
 def preconditioner_matrix_via_mapping_matrix_from(
     mapping_matrix: np.ndarray,
     regularization_matrix: np.ndarray,
diff --git a/autoarray/inversion/inversion/settings.py b/autoarray/inversion/inversion/settings.py
index 1f489e63a..8c83e5e95 100644
--- a/autoarray/inversion/inversion/settings.py
+++ b/autoarray/inversion/inversion/settings.py
@@ -12,8 +12,8 @@ def __init__(
         use_w_tilde: bool = True,
         use_positive_only_solver: bool = False,
         positive_only_maxiter: int = 5000,
+        force_edge_pixels_to_zeros: bool = False,
         no_regularization_add_to_curvature_diag: bool = True,
-        check_solution: Optional[bool] = None,
         use_w_tilde_numpy: bool = False,
         use_source_loop: bool = False,
         use_linear_operators: bool = False,
@@ -40,10 +40,6 @@ def __init__(
         no_regularization_add_to_curvature_diag
             When True, if a linear object in the inversion has no regularization, values of 1.0e-8 are added to the
             diagonal of its `curvature_matrix` to stablelize the linear algebra solver.
-        check_solution
-            If True, the `reconstruction` of the inversion is checked to ensure that no two source pixels have
-            numerically identical values, which indicates a spurious solution where the linear algebra solver
-            reconstructs every value as an identical value (for reasons that are currently not understood).
         use_w_tilde_numpy
             If True, the curvature_matrix is computed via numpy matrix multiplication (as opposed to numba functions
             which exploit sparsity to do the calculation normally in a more efficient way).
@@ -64,6 +60,7 @@ def __init__(
         self.use_positive_only_solver = use_positive_only_solver
         self.positive_only_maxiter = positive_only_maxiter
         self.use_linear_operators = use_linear_operators
+        self.force_edge_pixels_to_zeros = force_edge_pixels_to_zeros
         self.no_regularization_add_to_curvature_diag = (
             no_regularization_add_to_curvature_diag
         )
@@ -71,14 +68,3 @@ def __init__(
         self.maxiter = maxiter
         self.use_w_tilde_numpy = use_w_tilde_numpy
         self.use_source_loop = use_source_loop
-
-        self._check_solution = check_solution
-
-    @property
-    def check_solution(self):
-
-        if self._check_solution is None:
-
-            return conf.instance["general"]["inversion"]["check_solution_default"]
-
-        return self._check_solution
diff --git a/autoarray/inversion/mock/mock_mapper.py b/autoarray/inversion/mock/mock_mapper.py
index db44eba0a..9e9c69e60 100644
--- a/autoarray/inversion/mock/mock_mapper.py
+++ b/autoarray/inversion/mock/mock_mapper.py
@@ -11,6 +11,7 @@ def __init__(
         source_plane_data_grid=None,
         source_plane_mesh_grid=None,
         hyper_data=None,
+        edge_pixel_list=None,
         regularization=None,
         pix_sub_weights=None,
         mapping_matrix=None,
@@ -27,6 +28,8 @@ def __init__(
 
         super().__init__(mapper_grids=mapper_grids, regularization=regularization)
 
+        self._edge_pixel_list = edge_pixel_list
+
         self._pix_sub_weights = pix_sub_weights
 
         self._mapping_matrix = mapping_matrix
@@ -48,6 +51,10 @@ def params(self):
             return super().params
         return self._parameters
 
+    @property
+    def edge_pixel_list(self):
+        return self._edge_pixel_list
+
     @property
     def pix_sub_weights(self):
         return self._pix_sub_weights
diff --git a/autoarray/inversion/pixelization/mappers/abstract.py b/autoarray/inversion/pixelization/mappers/abstract.py
index 0dfd7144e..fa262c788 100644
--- a/autoarray/inversion/pixelization/mappers/abstract.py
+++ b/autoarray/inversion/pixelization/mappers/abstract.py
@@ -104,6 +104,10 @@ def source_plane_mesh_grid(self) -> Abstract2DMesh:
     def image_plane_mesh_grid(self) -> Grid2D:
         return self.mapper_grids.image_plane_mesh_grid
 
+    @property
+    def edge_pixel_list(self) -> List[int]:
+        return self.source_plane_mesh_grid.edge_pixel_list
+
     @property
     def hyper_data(self) -> np.ndarray:
         return self.mapper_grids.hyper_data
diff --git a/autoarray/inversion/pixelization/mappers/mapper_util.py b/autoarray/inversion/pixelization/mappers/mapper_util.py
index c09b485aa..bea6cb619 100644
--- a/autoarray/inversion/pixelization/mappers/mapper_util.py
+++ b/autoarray/inversion/pixelization/mappers/mapper_util.py
@@ -152,7 +152,7 @@ def pix_indexes_for_sub_slim_index_delaunay_from(
     delaunay_points,
 ) -> Tuple[np.ndarray, np.ndarray]:
     """
-    The indexes mappings between the sub pixels and Voronoi pixelization pixels.
+    The indexes mappings between the sub pixels and Voronoi mesh pixels.
     For Delaunay tessellation, most sub pixels should have contribution of 3 pixelization pixels. However,
     for those ones not belonging to any triangle, we link its value to its closest point.
 
diff --git a/autoarray/inversion/pixelization/mappers/voronoi.py b/autoarray/inversion/pixelization/mappers/voronoi.py
index b9e732d66..20bcc3519 100644
--- a/autoarray/inversion/pixelization/mappers/voronoi.py
+++ b/autoarray/inversion/pixelization/mappers/voronoi.py
@@ -136,10 +136,10 @@ def pix_sub_weights(self) -> PixSubWeights:
         neighbor of the Voronoi mesh's third (index 2) pixel.
 
         - `pix_indexes_for_sub_slim_index[0, 1] = 5`: The data's first (index 0) sub-pixel also maps to the natural
-        neighbor of the Voronoi pixelization's sixth (index 5) pixel.
+        neighbor of the Voronoi mesh's sixth (index 5) pixel.
 
         - `pix_indexes_for_sub_slim_index[0, 2] = 8`: The data's first (index 0) sub-pixel also maps to the natural
-        neighbor of the Voronoi pixelization's ninth (index 8) pixel.
+        neighbor of the Voronoi mesh's ninth (index 8) pixel.
 
         The interpolation weights of these multiple mappings are stored in the array `pix_weights_for_sub_slim_index`.
         """
@@ -222,7 +222,7 @@ def pix_sub_weights(self) -> PixSubWeights:
         where each data pixel maps to 3 Delaunay triangles with interpolation weights). The weights of multiple mappings
         are stored in the array `pix_weights_for_sub_slim_index`.
 
-        For this Voronoi pixelization each data sub-pixel maps to a single pixelization pixel, thus the second
+        For this Voronoi mesh each data sub-pixel maps to a single pixelization pixel, thus the second
         dimension of the array `pix_indexes_for_sub_slim_index` 1 and all entries in `pix_weights_for_sub_slim_index`
         are equal to 1.0.
 
diff --git a/autoarray/inversion/pixelization/mesh/mesh_util.py b/autoarray/inversion/pixelization/mesh/mesh_util.py
index 80b92c1b0..92c31ead7 100644
--- a/autoarray/inversion/pixelization/mesh/mesh_util.py
+++ b/autoarray/inversion/pixelization/mesh/mesh_util.py
@@ -295,6 +295,32 @@ def rectangular_central_neighbors(
     return neighbors, neighbors_sizes
 
 
+def rectangular_edge_pixel_list_from(neighbors: np.ndarray) -> List:
+    """
+    Returns a list of the 1D indices of all pixels on the edge of a rectangular pixelization.
+
+    This is computed by searching the `neighbors` array for pixels that have a neighbor with index -1, meaning there
+    is at least one neighbor from the 4 expected missing.
+
+    Parameters
+    ----------
+    neighbors
+        An array of dimensions [total_pixels, 4] which provides the index of all neighbors of every pixel in the
+        rectangular pixelization (entries of -1 correspond to no neighbor).
+
+    Returns
+    -------
+    A list of the 1D indices of all pixels on the edge of a rectangular pixelization.
+    """
+    edge_pixel_list = []
+
+    for i, neighbors in enumerate(neighbors):
+        if -1 in neighbors:
+            edge_pixel_list.append(i)
+
+    return edge_pixel_list
+
+
 @numba_util.jit()
 def delaunay_triangle_area_from(
     corner_0: Tuple[float, float],
@@ -420,7 +446,7 @@ def voronoi_neighbors_from(
     pixels: int, ridge_points: np.ndarray
 ) -> Tuple[np.ndarray, np.ndarray]:
     """
-    Returns the adjacent neighbors of every pixel on a Voronoi pixelization as an ndarray of shape
+    Returns the adjacent neighbors of every pixel on a Voronoi mesh as an ndarray of shape
     [total_pixels, voronoi_pixel_with_max_neighbors], using the `ridge_points` output from the `scipy.spatial.Voronoi()`
     object.
 
@@ -428,7 +454,7 @@ def voronoi_neighbors_from(
     an ndarray with the number of neighbors of every pixel, `neighbors_sizes`, which is iterated over when using
     the `neighbors` ndarray.
 
-    Indexing is defined in an arbritrary manner due to the irregular nature of a Voronoi pixelization.
+    Indexing is defined in an arbritrary manner due to the irregular nature of a Voronoi mesh.
 
     For example, if `neighbors[0,:] = [1, 5, 36, 2, -1, -1]`, this informs us that the first Voronoi pixel has
     4 neighbors which have indexes 1, 5, 36, 2. Correspondingly `neighbors_sizes[0] = 4`.
@@ -436,7 +462,7 @@ def voronoi_neighbors_from(
     Parameters
     ----------
     pixels
-        The number of pixels on the Voronoi pixelization.
+        The number of pixels on the Voronoi mesh.
     ridge_points
         Contains the information on every Voronoi source pixel and its neighbors.
 
@@ -444,7 +470,22 @@ def voronoi_neighbors_from(
     -------
     The arrays containing the 1D index of every pixel's neighbors and the number of neighbors that each pixel has.
     """
-
+    """
+    Returns the neighbors and total number of neighbors of every cell on a Voronoi mesh. 
+    
+    Neighbors are returned as an ndarray of shape [total_pixels, max_neighbors_in_a_given_voronoi pixel], where 
+    entries have values of -1 if the pixel has no neighbor.
+    
+    The number of neighbors of every pixel is also returned as an ndarray of shape [total_pixels], where the values
+    are integers between 0 and the total neighbors in a given Voronoi pixel.
+    
+    Parameters
+    ----------
+    pixels
+        The number of pixels on the Voronoi mesh
+    ridge_points
+        Contains the information on every Voronoi pixel and its neighbors.
+    """
     neighbors_sizes = np.zeros(shape=(pixels))
 
     for ridge_index in range(ridge_points.shape[0]):
@@ -467,11 +508,33 @@ def voronoi_neighbors_from(
     return neighbors, neighbors_sizes
 
 
+# @numba_util.jit()
+def voronoi_edge_pixels_from(regions: np.ndarray, point_region: np.ndarray) -> List:
+    """
+    Returns the edge pixels of a Voronoi mesh, where the edge pixels are defined as those pixels which are on the
+    edge of the Voronoi diagram.
+
+    Parameters
+    ----------
+    regions
+        Indices of the Voronoi vertices forming each Voronoi region, where -1 indicates vertex outside the Voronoi
+        diagram.
+    """
+
+    voronoi_edge_pixel_list = []
+
+    for index, i in enumerate(point_region):
+        if -1 in regions[i]:
+            voronoi_edge_pixel_list.append(index)
+
+    return voronoi_edge_pixel_list
+
+
 def voronoi_revised_from(
     voronoi: scipy.spatial.Voronoi,
 ) -> Union[List[Tuple], np.ndarray]:
     """
-    To plot a Voronoi pixelization using the `matplotlib.fill()` function a revised Voronoi pixelization must be
+    To plot a Voronoi mesh using the `matplotlib.fill()` function a revised Voronoi mesh must be
      computed, where 2D infinite voronoi regions are converted to finite 2D regions.
 
     This function returns a list of tuples containing the indices of the vertices of each revised Voronoi cell and
diff --git a/autoarray/inversion/pixelization/mesh/voronoi.py b/autoarray/inversion/pixelization/mesh/voronoi.py
index 79067a84b..5f5b5a693 100644
--- a/autoarray/inversion/pixelization/mesh/voronoi.py
+++ b/autoarray/inversion/pixelization/mesh/voronoi.py
@@ -52,7 +52,7 @@ def mesh_grid_from(
     ) -> Mesh2DVoronoi:
         """
         Return the Voronoi `source_plane_mesh_grid` as a `Mesh2DVoronoi` object, which provides additional
-        functionality for performing operations that exploit the geometry of a Voronoi pixelization.
+        functionality for performing operations that exploit the geometry of a Voronoi mesh.
 
         The array `sparse_index_for_slim_index` encodes the closest source pixel of every pixel on the
         (full resolution) sub image-plane grid. This is used for efficiently pairing every image-plane pixel to its
@@ -232,7 +232,7 @@ def image_plane_mesh_grid_from(
         masked 2D data (see `Grid2DSparse.from_grid_and_unmasked_2d_grid_shape()`).
 
         The `data_pixelization_grid` is transformed to the `source_plane_mesh_grid`, and it is these (y,x) values
-        which then act the centres of the Voronoi pixelization's pixels.
+        which then act the centres of the Voronoi mesh's pixels.
 
         For a `VoronoiBrightnessImage` this grid is computed by applying a KMeans clustering algorithm to the masked
         data's values, where these values are reweighted by the `hyper_data` so that the algorithm can adapt to
diff --git a/autoarray/inversion/regularization/abstract.py b/autoarray/inversion/regularization/abstract.py
index cd2bcde23..137650d0d 100644
--- a/autoarray/inversion/regularization/abstract.py
+++ b/autoarray/inversion/regularization/abstract.py
@@ -87,7 +87,7 @@ def __init__(self):
         Whilst the example above used a square-grid with regularization to the right and downwards, this matrix \
         formalism can be extended to describe regularization in more directions (e.g. upwards, to the left).
 
-        It can also describe irpixelizations, e.g. an irVoronoi pixelization, where a B matrix is \
+        It can also describe irpixelizations, e.g. an irVoronoi mesh, where a B matrix is \
         computed for every shared Voronoi vertex of each Voronoi pixel. The number of B matrices is now equal to the \
         number of Voronoi vertices in the pixel with the most Voronoi vertices. However, we describe below a scheme to \
         compute this solution more efficiently.
diff --git a/autoarray/structures/mesh/delaunay_2d.py b/autoarray/structures/mesh/delaunay_2d.py
index bdd90b24a..a6b7f9012 100644
--- a/autoarray/structures/mesh/delaunay_2d.py
+++ b/autoarray/structures/mesh/delaunay_2d.py
@@ -1,5 +1,5 @@
 import numpy as np
-from typing import Optional, Tuple
+from typing import List, Optional, Tuple
 
 from autoconf import cached_property
 
@@ -18,7 +18,7 @@ def neighbors(self) -> Neighbors:
 
         see `Neighbors` for a complete description of the neighboring scheme.
 
-        The neighbors of a Voronoi pixelization are computed using the `ridge_points` attribute of the scipy `Voronoi`
+        The neighbors of a Voronoi mesh are computed using the `ridge_points` attribute of the scipy `Voronoi`
         object, as described in the method `mesh_util.voronoi_neighbors_from`.
         """
         indptr, indices = self.delaunay.vertex_neighbor_vertices
diff --git a/autoarray/structures/mesh/rectangular_2d.py b/autoarray/structures/mesh/rectangular_2d.py
index 5641e515e..3702373b6 100644
--- a/autoarray/structures/mesh/rectangular_2d.py
+++ b/autoarray/structures/mesh/rectangular_2d.py
@@ -1,6 +1,6 @@
 import numpy as np
 from scipy.interpolate import griddata
-from typing import Optional, Tuple
+from typing import List, Optional, Tuple
 
 from autoconf import cached_property
 
@@ -132,6 +132,10 @@ def neighbors(self) -> Neighbors:
 
         return Neighbors(arr=neighbors.astype("int"), sizes=sizes.astype("int"))
 
+    @cached_property
+    def edge_pixel_list(self) -> List:
+        return mesh_util.rectangular_edge_pixel_list_from(neighbors=self.neighbors)
+
     @property
     def pixels(self) -> int:
         """
diff --git a/autoarray/structures/mesh/triangulation_2d.py b/autoarray/structures/mesh/triangulation_2d.py
index dc40a0772..3c244ef0f 100644
--- a/autoarray/structures/mesh/triangulation_2d.py
+++ b/autoarray/structures/mesh/triangulation_2d.py
@@ -9,6 +9,7 @@
 from autoarray.structures.mesh.abstract_2d import Abstract2DMesh
 
 from autoarray import exc
+from autoarray.inversion.pixelization.mesh import mesh_util
 from autoarray.structures.grids import grid_2d_util
 
 
@@ -39,7 +40,7 @@ def __new__(
         The input `grid` of source pixel centres is ordered arbitrarily, given that there is no regular pattern
         for a Delaunay triangulation and Voronoi mesh's indexing to follow.
 
-        This class is used in conjuction with the `inversion/pixelizations` package to create Voronoi pixelizations
+        This class is used in conjuction with the `inversion/pixelizations` package to create Voronoi meshs
         and mappers that perform an `Inversion`.
 
         Parameters
@@ -48,7 +49,7 @@ def __new__(
             The grid of (y,x) coordinates corresponding to the Delaunay triangle corners and Voronoi pixel centres.
         nearest_pixelization_index_for_slim_index
             When a Voronoi grid is used to create a mapper and inversion, there are mappings between the `data` pixels
-            and Voronoi pixelization. This array contains these mappings and it is used to speed up the creation of the
+            and Voronoi mesh. This array contains these mappings and it is used to speed up the creation of the
             mapper.
         """
 
@@ -139,6 +140,16 @@ def voronoi(self) -> scipy.spatial.Voronoi:
         except (ValueError, OverflowError, scipy.spatial.qhull.QhullError) as e:
             raise exc.MeshException() from e
 
+    @cached_property
+    def edge_pixel_list(self) -> List:
+        """
+        Returns a list of the Voronoi pixel indexes that are on the edge of the mesh.
+        """
+
+        return mesh_util.voronoi_edge_pixels_from(
+            regions=self.voronoi.regions, point_region=self.voronoi.point_region
+        )
+
     @cached_property
     def split_cross(self) -> np.ndarray:
         """
@@ -227,6 +238,6 @@ def origin(self) -> Tuple[float, float]:
     @property
     def pixels(self) -> int:
         """
-        The total number of pixels in the Voronoi pixelization.
+        The total number of pixels in the Voronoi mesh.
         """
         return self.shape[0]
diff --git a/autoarray/structures/mesh/voronoi_2d.py b/autoarray/structures/mesh/voronoi_2d.py
index 56298e76a..f46b03d6d 100644
--- a/autoarray/structures/mesh/voronoi_2d.py
+++ b/autoarray/structures/mesh/voronoi_2d.py
@@ -1,6 +1,6 @@
 import numpy as np
 from scipy.interpolate import griddata
-from typing import Optional, Tuple
+from typing import List, Optional, Tuple
 
 from autoconf import cached_property
 
@@ -20,9 +20,10 @@ def neighbors(self) -> Neighbors:
 
         see `Neighbors` for a complete description of the neighboring scheme.
 
-        The neighbors of a Voronoi pixelization are computed using the `ridge_points` attribute of the scipy `Voronoi`
+        The neighbors of a Voronoi mesh are computed using the `ridge_points` attribute of the scipy `Voronoi`
         object, as described in the method `mesh_util.voronoi_neighbors_from`.
         """
+
         neighbors, sizes = mesh_util.voronoi_neighbors_from(
             pixels=self.pixels, ridge_points=np.asarray(self.voronoi.ridge_points)
         )
diff --git a/autoarray/util/nn/README.md b/autoarray/util/nn/README.md
index 83fc1d083..7652af031 100644
--- a/autoarray/util/nn/README.md
+++ b/autoarray/util/nn/README.md
@@ -4,7 +4,7 @@ Natural Neighbour
 (Natural Neighbours interpolation for PyAutoLens)
 
 If you want to use the `VoronoiNN` pixelization, which applies natural neighbor 
-interpolation (https://en.wikipedia.org/wiki/Natural_neighbor_interpolation) to a Voronoi pixelization you must 
+interpolation (https://en.wikipedia.org/wiki/Natural_neighbor_interpolation) to a Voronoi mesh you must 
 install this C package. 
 
 This currently requires that PyAutoLens is built from source, e.g. via cloning PyAutoLens and its parent packagees 
diff --git a/test_autoarray/inversion/inversion/test_abstract.py b/test_autoarray/inversion/inversion/test_abstract.py
index e3d77bc3e..b2bc7fb09 100644
--- a/test_autoarray/inversion/inversion/test_abstract.py
+++ b/test_autoarray/inversion/inversion/test_abstract.py
@@ -57,6 +57,20 @@ def test__index_range_list_from():
     assert inversion.param_range_list_from(cls=aa.AbstractMapper) == [[2, 3]]
 
 
+def test__mapper_edge_pixel_list():
+
+    inversion = aa.m.MockInversion(
+        linear_obj_list=[
+            aa.m.MockLinearObj(parameters=3, regularization=None),
+            aa.m.MockMapper(parameters=4, edge_pixel_list=[0, 2], regularization=None),
+            aa.m.MockLinearObj(parameters=7, regularization=None),
+            aa.m.MockMapper(parameters=4, edge_pixel_list=[0, 2], regularization=None),
+        ]
+    )
+
+    assert inversion.mapper_edge_pixel_list == [3, 5, 14, 16]
+
+
 def test__no_regularization_index_list():
 
     inversion = aa.m.MockInversion(
@@ -141,13 +155,13 @@ def test__curvature_matrix__via_w_tilde__identical_to_mapping():
     inversion_w_tilde = aa.Inversion(
         dataset=masked_imaging,
         linear_obj_list=[mapper_0, mapper_1],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     inversion_mapping = aa.Inversion(
         dataset=masked_imaging,
         linear_obj_list=[mapper_0, mapper_1],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert inversion_w_tilde.curvature_matrix == pytest.approx(
@@ -212,13 +226,13 @@ def test__curvature_matrix_via_w_tilde__includes_source_interpolation__identical
     inversion_w_tilde = aa.Inversion(
         dataset=masked_imaging,
         linear_obj_list=[mapper_0, mapper_1],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     inversion_mapping = aa.Inversion(
         dataset=masked_imaging,
         linear_obj_list=[mapper_0, mapper_1],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert inversion_w_tilde.curvature_matrix == pytest.approx(
@@ -244,6 +258,33 @@ def test__curvature_reg_matrix_reduced():
     ).all()
 
 
+def test__curvature_reg_matrix_solver__edge_pixels_set_to_zero():
+
+    curvature_reg_matrix = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])
+
+    linear_obj_list = [
+        aa.m.MockMapper(parameters=3, regularization=None, edge_pixel_list=[0])
+    ]
+
+    inversion = aa.m.MockInversion(
+        linear_obj_list=linear_obj_list,
+        curvature_reg_matrix=curvature_reg_matrix,
+        settings=aa.SettingsInversion(force_edge_pixels_to_zeros=True),
+    )
+
+    curvature_reg_matrix = np.array(
+        [
+            [0.0, 2.0, 3.0],
+            [0.0, 5.0, 6.0],
+            [0.0, 8.0, 9.0],
+        ]
+    )
+
+    assert inversion.curvature_reg_matrix_solver == pytest.approx(
+        curvature_reg_matrix, 1.0e-4
+    )
+
+
 def test__regularization_matrix():
 
     reg_0 = aa.m.MockRegularization(regularization_matrix=np.ones((2, 2)))
diff --git a/test_autoarray/inversion/inversion/test_factory.py b/test_autoarray/inversion/inversion/test_factory.py
index ab0a9f764..f554c1ece 100644
--- a/test_autoarray/inversion/inversion/test_factory.py
+++ b/test_autoarray/inversion/inversion/test_factory.py
@@ -18,7 +18,7 @@ def test__inversion_imaging__via_linear_obj_func_list(masked_imaging_7x7_no_blur
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[linear_obj],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.m.MockLinearObjFuncList)
@@ -31,7 +31,7 @@ def test__inversion_imaging__via_linear_obj_func_list(masked_imaging_7x7_no_blur
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[linear_obj],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.m.MockLinearObjFuncList)
@@ -48,7 +48,7 @@ def test__inversion_imaging__via_linear_obj_func_list(masked_imaging_7x7_no_blur
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[linear_obj],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.m.MockLinearObjFuncList)
@@ -67,7 +67,7 @@ def test__inversion_imaging__via_mapper(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[rectangular_mapper_7x7_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangularNoInterp)
@@ -78,7 +78,7 @@ def test__inversion_imaging__via_mapper(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[rectangular_mapper_7x7_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangularNoInterp)
@@ -89,7 +89,7 @@ def test__inversion_imaging__via_mapper(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[delaunay_mapper_9_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperDelaunay)
@@ -100,7 +100,7 @@ def test__inversion_imaging__via_mapper(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[delaunay_mapper_9_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperDelaunay)
@@ -111,7 +111,7 @@ def test__inversion_imaging__via_mapper(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[voronoi_mapper_9_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperVoronoiNoInterp)
@@ -122,7 +122,7 @@ def test__inversion_imaging__via_mapper(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[voronoi_mapper_9_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperVoronoiNoInterp)
@@ -148,7 +148,7 @@ def test__inversion_imaging__via_regularizations(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[mapper],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperDelaunay)
@@ -163,7 +163,7 @@ def test__inversion_imaging__via_regularizations(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[mapper],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperDelaunay)
@@ -178,7 +178,7 @@ def test__inversion_imaging__via_regularizations(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[mapper],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperVoronoiNoInterp)
@@ -191,7 +191,7 @@ def test__inversion_imaging__via_regularizations(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[mapper],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperVoronoiNoInterp)
@@ -210,7 +210,7 @@ def test__inversion_imaging__via_regularizations(
         inversion = aa.Inversion(
             dataset=masked_imaging_7x7_no_blur,
             linear_obj_list=[mapper],
-            settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+            settings=aa.SettingsInversion(use_w_tilde=True),
         )
 
         assert isinstance(inversion.linear_obj_list[0], aa.MapperVoronoiNoInterp)
@@ -225,7 +225,7 @@ def test__inversion_imaging__via_regularizations(
         inversion = aa.Inversion(
             dataset=masked_imaging_7x7_no_blur,
             linear_obj_list=[mapper],
-            settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+            settings=aa.SettingsInversion(use_w_tilde=True),
         )
 
         assert isinstance(inversion.linear_obj_list[0], aa.MapperVoronoi)
@@ -240,7 +240,7 @@ def test__inversion_imaging__via_regularizations(
         inversion = aa.Inversion(
             dataset=masked_imaging_7x7_no_blur,
             linear_obj_list=[mapper],
-            settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+            settings=aa.SettingsInversion(use_w_tilde=True),
         )
 
         assert isinstance(inversion.linear_obj_list[0], aa.MapperVoronoi)
@@ -276,7 +276,6 @@ def test__inversion_imaging__via_linear_obj_func_and_mapper(
         linear_obj_list=[linear_obj, rectangular_mapper_7x7_3x3],
         settings=aa.SettingsInversion(
             use_w_tilde=False,
-            check_solution=False,
             no_regularization_add_to_curvature_diag=False,
         ),
     )
@@ -294,6 +293,57 @@ def test__inversion_imaging__via_linear_obj_func_and_mapper(
     assert inversion.mapped_reconstructed_image == pytest.approx(np.ones(9), 1.0e-4)
 
 
+def test__inversion_imaging__via_linear_obj_func_and_mapper__force_edge_pixels_to_zero(
+    masked_imaging_7x7_no_blur,
+    voronoi_mapper_9_3x3,
+):
+
+    mask = masked_imaging_7x7_no_blur.mask
+
+    grid = aa.Grid2D.from_mask(mask=mask)
+
+    linear_obj = aa.m.MockLinearObj(
+        parameters=1,
+        grid=grid,
+        mapping_matrix=np.full(fill_value=0.5, shape=(9, 1)),
+        regularization=None,
+    )
+
+    inversion = aa.Inversion(
+        dataset=masked_imaging_7x7_no_blur,
+        linear_obj_list=[linear_obj, voronoi_mapper_9_3x3],
+        settings=aa.SettingsInversion(
+            use_w_tilde=False,
+            no_regularization_add_to_curvature_diag=False,
+            force_edge_pixels_to_zeros=True,
+        ),
+    )
+
+    assert isinstance(inversion.linear_obj_list[0], aa.m.MockLinearObj)
+    assert isinstance(inversion.linear_obj_list[1], aa.MapperVoronoiNoInterp)
+    assert isinstance(inversion, aa.InversionImagingMapping)
+    assert (inversion.curvature_reg_matrix_solver[:, 1] == np.zeros(shape=(10,))).all()
+
+    inversion = aa.Inversion(
+        dataset=masked_imaging_7x7_no_blur,
+        linear_obj_list=[linear_obj, voronoi_mapper_9_3x3],
+        settings=aa.SettingsInversion(
+            use_w_tilde=False,
+            use_positive_only_solver=True,
+            no_regularization_add_to_curvature_diag=False,
+            force_edge_pixels_to_zeros=True,
+        ),
+    )
+
+    assert isinstance(inversion.linear_obj_list[0], aa.m.MockLinearObj)
+    assert isinstance(inversion.linear_obj_list[1], aa.MapperVoronoiNoInterp)
+    assert isinstance(inversion, aa.InversionImagingMapping)
+    assert (inversion.curvature_reg_matrix_solver[:, 1] == np.zeros(shape=(10,))).all()
+    assert inversion.reconstruction == pytest.approx(
+        np.array([2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), 1.0e-4
+    )
+
+
 def test__inversion_imaging__compare_mapping_and_w_tilde_values(
     masked_imaging_7x7, voronoi_mapper_9_3x3
 ):
@@ -346,7 +396,7 @@ def test__inversion_imaging__linear_obj_func_and_non_func_give_same_terms(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[linear_obj, rectangular_mapper_7x7_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     masked_imaging_7x7_no_blur = copy.copy(masked_imaging_7x7_no_blur)
@@ -358,7 +408,7 @@ def test__inversion_imaging__linear_obj_func_and_non_func_give_same_terms(
     inversion_no_linear_func = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[rectangular_mapper_7x7_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert inversion.regularization_term == pytest.approx(
@@ -400,13 +450,13 @@ def test__inversion_imaging__linear_obj_func_with_w_tilde(
     inversion_mapping = aa.Inversion(
         dataset=masked_imaging_7x7,
         linear_obj_list=[linear_obj, rectangular_mapper_7x7_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     inversion_w_tilde = aa.Inversion(
         dataset=masked_imaging_7x7,
         linear_obj_list=[linear_obj, rectangular_mapper_7x7_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     assert inversion_mapping.data_vector == pytest.approx(
@@ -437,7 +487,7 @@ def test__inversion_imaging__linear_obj_func_with_w_tilde(
             voronoi_mapper_9_3x3,
             linear_obj_2,
         ],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     inversion_w_tilde = aa.Inversion(
@@ -450,7 +500,7 @@ def test__inversion_imaging__linear_obj_func_with_w_tilde(
             voronoi_mapper_9_3x3,
             linear_obj_2,
         ],
-        settings=aa.SettingsInversion(use_w_tilde=True, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=True),
     )
 
     assert inversion_mapping.data_vector == pytest.approx(
@@ -471,7 +521,7 @@ def test__inversion_interferometer__via_mapper(
     inversion = aa.Inversion(
         dataset=interferometer_7_no_fft,
         linear_obj_list=[rectangular_mapper_7x7_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperRectangularNoInterp)
@@ -486,7 +536,7 @@ def test__inversion_interferometer__via_mapper(
     inversion = aa.Inversion(
         dataset=interferometer_7_no_fft,
         linear_obj_list=[delaunay_mapper_9_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperDelaunay)
@@ -503,7 +553,7 @@ def test__inversion_interferometer__via_mapper(
     inversion = aa.Inversion(
         dataset=interferometer_7_no_fft,
         linear_obj_list=[voronoi_mapper_9_3x3],
-        settings=aa.SettingsInversion(use_w_tilde=False, check_solution=False),
+        settings=aa.SettingsInversion(use_w_tilde=False),
     )
 
     assert isinstance(inversion.linear_obj_list[0], aa.MapperVoronoiNoInterp)
@@ -526,7 +576,6 @@ def test__inversion_matrices__x2_mappers(
     inversion = aa.Inversion(
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[rectangular_mapper_7x7_3x3, voronoi_mapper_9_3x3],
-        settings=aa.SettingsInversion(check_solution=False),
     )
 
     assert (
@@ -604,7 +653,7 @@ def test__inversion_imaging__positive_only_solver(masked_imaging_7x7_no_blur):
         dataset=masked_imaging_7x7_no_blur,
         linear_obj_list=[linear_obj],
         settings=aa.SettingsInversion(
-            use_w_tilde=False, check_solution=False, use_positive_only_solver=True
+            use_w_tilde=False, use_positive_only_solver=True
         ),
     )
 
diff --git a/test_autoarray/structures/mesh/test_delaunay.py b/test_autoarray/structures/mesh/test_delaunay.py
index 214aeedc7..96956e923 100644
--- a/test_autoarray/structures/mesh/test_delaunay.py
+++ b/test_autoarray/structures/mesh/test_delaunay.py
@@ -43,6 +43,32 @@ def test__mesh_areas():
     )
 
 
+def test__edge_pixel_list():
+
+    grid = np.array(
+        [
+            [1.0, -1.0],
+            [1.0, 0.0],
+            [1.0, 1.0],
+            [0.0, -1.0],
+            [0.0, 0.0],
+            [0.0, 1.0],
+            [-1.0, -1.0],
+            [-1.0, 0.0],
+            [-1.0, 1.0],
+        ]
+    )
+
+    mesh = aa.Mesh2DDelaunay(values=grid)
+
+    print(mesh.delaunay.neighbors)
+    print(mesh.delaunay.vertex_neighbor_vertices)
+    print(mesh.neighbors)
+    print(mesh.edge_pixel_list)
+
+    assert mesh.edge_pixel_list == [0, 1, 2, 3, 5, 6, 7, 8]
+
+
 def test__interpolated_array_from():
 
     grid = aa.Grid2D.no_mask(
diff --git a/test_autoarray/structures/mesh/test_rectangular.py b/test_autoarray/structures/mesh/test_rectangular.py
index 45925fd4f..b1fe78bd5 100644
--- a/test_autoarray/structures/mesh/test_rectangular.py
+++ b/test_autoarray/structures/mesh/test_rectangular.py
@@ -24,6 +24,28 @@ def test__neighbors__compare_to_mesh_util():
     assert (mesh.neighbors.sizes == neighbors_sizes_util).all()
 
 
+def test__edge_pixel_list():
+    grid = np.array(
+        [
+            [-1.0, -1.0],
+            [-1.0, 0.0],
+            [-1.0, 1.0],
+            [0.0, -1.0],
+            [0.0, 0.0],
+            [0.0, 1.0],
+            [1.0, -1.0],
+            [1.0, 0.0],
+            [1.0, 1.0],
+        ]
+    )
+
+    mesh = aa.Mesh2DRectangular.overlay_grid(
+        shape_native=(3, 3), grid=grid, buffer=1e-8
+    )
+
+    assert mesh.edge_pixel_list == [0, 1, 2, 3, 5, 6, 7, 8]
+
+
 def test__shape_native_and_pixel_scales():
     grid = np.array(
         [
diff --git a/test_autoarray/structures/mesh/test_voronoi.py b/test_autoarray/structures/mesh/test_voronoi.py
index 6ce47d617..00c9b25f9 100644
--- a/test_autoarray/structures/mesh/test_voronoi.py
+++ b/test_autoarray/structures/mesh/test_voronoi.py
@@ -88,6 +88,27 @@ def test__mesh_grid__attributes():
     ).all()
 
 
+def test__edge_pixel_list():
+
+    grid = np.array(
+        [
+            [1.0, -1.0],
+            [1.0, 0.0],
+            [1.0, 1.0],
+            [0.0, -1.0],
+            [0.0, 0.0],
+            [0.0, 1.0],
+            [-1.0, -1.0],
+            [-1.0, 0.0],
+            [-1.0, 1.0],
+        ]
+    )
+
+    mesh = aa.Mesh2DVoronoi(values=grid)
+
+    assert mesh.edge_pixel_list == [0, 1, 2, 3, 5, 6, 7, 8]
+
+
 def test__from_unmasked_sparse_shape_and_grid():
 
     mask = aa.Mask2D(