@@ -150,6 +150,25 @@ class InterpolatorKNearestNeighbor(InterpolatorDelaunay):
     # Excluding them would change pixel signals broadly and add a jump in the likelihood (PyAutoArray#627).
     exclude_zero_weight_mappings = False
 
+    @cached_property
+    def mesh_geometry(self):
+        """
+        The kNN geometry: Delaunay vertex geometry whose
+        `areas_for_magnification` integrates this interpolator's own basis
+        functions (see `MeshGeometryKNearestNeighbor`).
+        """
+        from autoarray.inversion.mesh.mesh_geometry.delaunay import (
+            MeshGeometryKNearestNeighbor,
+        )
+
+        return MeshGeometryKNearestNeighbor(
+            mesh=self.mesh,
+            mesh_grid=self.mesh_grid,
+            data_grid=self.data_grid,
+            interpolator_cls=type(self),
+            xp=self._xp,
+        )
+
     @property
     def dual_areas(self):
         """
 
@@ -271,3 +271,87 @@ def areas_for_magnification(self) -> np.ndarray:
             scipy.spatial.Delaunay(mesh_grid_xy).simplices,
             xp=np,
         )
+
+
+# Quadrature cells per axis for `MeshGeometryKNearestNeighbor.areas_for_magnification`.
+# On a 144-vertex mesh the identity-lens magnification moved by < 3e-5 between
+# 200, 400 and 800 cells per axis (PyAutoArray#665); 256 costs ~1 s there.
+KNN_AREA_QUADRATURE_PIXELS: int = 256
+
+
+class MeshGeometryKNearestNeighbor(MeshGeometryDelaunay):
+    def __init__(self, mesh, mesh_grid, data_grid, interpolator_cls, xp=np, **kwargs):
+        """
+        The geometry of the kNN meshes (`KNearestNeighbor`, `KNNBarycentric`):
+        Delaunay vertex geometry, but magnification areas that integrate the
+        kNN interpolant itself.
+
+        Parameters
+        ----------
+        interpolator_cls
+            The kNN interpolator class that built this geometry; its basis
+            functions are evaluated on a quadrature grid to compute
+            `areas_for_magnification`.
+        """
+        super().__init__(
+            mesh=mesh,
+            mesh_grid=mesh_grid,
+            data_grid=data_grid,
+            dual_areas=None,
+            xp=xp,
+            **kwargs,
+        )
+
+        self.interpolator_cls = interpolator_cls
+
+    @cached_property
+    def areas_for_magnification(self) -> np.ndarray:
+        """
+        The integral of every vertex's interpolation basis function over the
+        convex hull of the mesh, so `sum(reconstruction * areas_for_magnification)`
+        is the integral of the reconstructed source.
+
+        The Delaunay barycentric dual areas are these integrals only for the linear
+        barycentric interpolant. The Wendland-kNN and 3-nearest barycentric
+        interpolants are not linear on the Delaunay triangles, and weighting by the
+        dual areas biased their magnification by up to ~6% (PyAutoArray#665).
+
+        The integral is a midpoint quadrature: the interpolator's own mapping
+        weights are evaluated at the centres of a uniform
+        `KNN_AREA_QUADRATURE_PIXELS`^2 grid over the vertices' bounding box, the
+        centres outside the convex hull are dropped, and the weights times the
+        cell area are summed per vertex. The total is rescaled to the hull area,
+        which removes the cells straddling the hull edge from the sum (the same
+        domain the Delaunay dual areas tile).
+
+        Computed host-side on first access, like the Delaunay dual-area fallback.
+        """
+        import scipy.spatial
+
+        nodes = np.asarray(self.mesh_grid_xy)
+        lo = nodes.min(axis=0)
+        hi = nodes.max(axis=0)
+
+        n = KNN_AREA_QUADRATURE_PIXELS
+        u = (np.arange(n) + 0.5) / n
+        grid_0, grid_1 = np.meshgrid(
+            lo[0] + u * (hi[0] - lo[0]), lo[1] + u * (hi[1] - lo[1]), indexing="ij"
+        )
+        points = np.stack([grid_0.ravel(), grid_1.ravel()], axis=1)
+        points = points[scipy.spatial.Delaunay(nodes).find_simplex(points) >= 0]
+
+        interpolator = self.interpolator_cls(
+            mesh=self.mesh,
+            mesh_grid=self.mesh_grid,
+            data_grid=points,
+            xp=np,
+        )
+        mappings, _, weights = interpolator._mappings_sizes_weights
+
+        areas = np.bincount(
+            np.asarray(mappings).ravel(),
+            weights=np.asarray(weights).ravel(),
+            minlength=nodes.shape[0],
+        )
+
+        return areas * scipy.spatial.ConvexHull(nodes).volume / np.sum(areas)
@@ -76,11 +76,18 @@
 #     not the integral of the node's hat once the CDF bends across the cell (the
 #     rank CDF's staircase on tied coordinates; a steep adapt-image CDF).
 #
-# Documented approximations (flagged in PyAutoArray#659; phase 4b):
-#   - DelaunayNN / KNearestNeighbor / KNNBarycentric inherit the Delaunay dual
-#     areas, which are exact only for the linear interpolant; for the Sibson,
-#     Wendland-kNN and 3-nearest-barycentric interpolants they are an approximation
-#     [-4.3e-3, +1.6e-2, -4.4e-2].
+#   - KNearestNeighbor / KNNBarycentric: the integral of each vertex's kNN basis
+#     function over the convex hull, by midpoint quadrature of the interpolator's
+#     own weights (`MeshGeometryKNearestNeighbor`, PyAutoArray#665)
+#     [+1.3e-5, +1.8e-3; KNNBarycentric's clipped weights are discontinuous, so
+#     the coarse sub-pixel `Q` is the limit there]. The Delaunay dual areas they
+#     used before measured [+1.6e-2, -4.4e-2]. The quadrature integrates whatever
+#     neighbours the interpolator selects, so it stays exact with or without the
+#     open partial-block neighbour bug PyAutoArray#609.
+#
+# Documented approximation:
+#   - DelaunayNN inherits the Delaunay dual areas, exact only for the linear
+#     interpolant; for its Sibson interpolant they are an approximation [-4.3e-3].
 MAGNIFICATION_TOLERANCE = {
     "RectangularUniform": 2.5e-3,
     "RectangularBilinearAdaptDensity": 1.0e-3,
@@ -89,8 +96,8 @@
     "RectangularRTUAdaptImage": 1.0e-2,
     "Delaunay": 3.0e-4,
     "DelaunayNN": 1.0e-2,
-    "KNearestNeighbor": 3.0e-2,
-    "KNNBarycentric": 7.0e-2,
+    "KNearestNeighbor": 1.0e-4,
+    "KNNBarycentric": 4.0e-3,
 }
 
 RECTANGULAR_ADAPTIVE = [