diff --git a/yt/utilities/lib/image_utilities.pyx b/yt/utilities/lib/image_utilities.pyx
index 686ed7c48d..d15e0f09e7 100644
--- a/yt/utilities/lib/image_utilities.pyx
+++ b/yt/utilities/lib/image_utilities.pyx
@@ -11,12 +11,46 @@ import numpy as np
 
 cimport numpy as np
 cimport cython
-from libc.math cimport ceil, floor, log2, sqrt
+from libc.math cimport atan2, ceil, floor, log2, sqrt, M_PI, M_PI_2
 from libc.stdlib cimport free, malloc
 
 from yt.utilities.lib.fp_utils cimport iclip, imin, imax, fclip, fmin, fmax
 from cython.parallel import prange, parallel
 
+cdef np.float64_t twopi = 2 * M_PI
+
+@cython.boundscheck(False)
+cdef inline void _add_point_to_greyscale_image(
+    np.float64_t[:, ::1] buffer,
+    np.uint8_t[:, ::1] buffer_mask,
+    np.float64_t px,
+    np.float64_t py,
+    np.float64_t pv,
+    int xs,
+    int ys
+) noexcept nogil:
+    cdef int i, j
+    j = (<int> (xs * px)) % xs
+    i = (<int> (ys * py)) % ys
+    if (i < 0) or (i >= buffer.shape[0]) or (j < 0) or (j >= buffer.shape[1]):
+        # some particles might intersect the image buffer
+        # but actually be centered out of bounds. Skip those.
+        # see https://github.com/yt-project/yt/issues/4603
+        return
+    buffer[i, j] += pv
+    buffer_mask[i, j] = 1
+
+
+
+cdef inline np.float64_t  _wrap_dist(np.float64_t dx) nogil:
+    if dx > 0.5:
+        return 1 - dx
+    elif dx < -0.5:
+        return 1 + dx
+    else:
+        return dx
+
+
 @cython.wraparound(False)
 def add_points_to_greyscale_image(
         np.ndarray[np.float64_t, ndim=2] buffer,
@@ -24,21 +58,16 @@ def add_points_to_greyscale_image(
         np.ndarray[np.float64_t, ndim=1] px,
         np.ndarray[np.float64_t, ndim=1] py,
         np.ndarray[np.float64_t, ndim=1] pv):
-    cdef int i, j, pi
+    cdef int pi
     cdef int npx = px.shape[0]
     cdef int xs = buffer.shape[0]
     cdef int ys = buffer.shape[1]
-    for pi in range(npx):
-        j = <int> (xs * px[pi])
-        i = <int> (ys * py[pi])
-        if (i < 0) or (i >= buffer.shape[0]) or (j < 0) or (j >= buffer.shape[1]):
-            # some particles might intersect the image buffer
-            # but actually be centered out of bounds. Skip those.
-            # see https://github.com/yt-project/yt/issues/4603
-            continue
 
-        buffer[i, j] += pv[pi]
-        buffer_mask[i, j] = 1
+    cdef np.float64_t[:, ::1] buffer_view = buffer
+    cdef np.uint8_t[:, ::1] buffer_mask_view = buffer_mask
+
+    for pi in range(npx):
+        _add_point_to_greyscale_image(buffer_view, buffer_mask_view, px[pi], py[pi], pv[pi], xs, ys)
     return
 
 cdef inline int ij2idx(const int i, const int j, const int Nx) noexcept nogil:
@@ -461,3 +490,373 @@ def add_rgba_points_to_image(
         for k in range(4):
             buffer[i, j, k] += rgba[pi, k]
     return
+
+@cython.boundscheck(False)
+cdef np.ndarray[np.float64_t, ndim=2] sample_tetrahedron(int sample):
+    cdef int i
+    cdef np.float64_t a, b, c
+
+    cdef np.ndarray[np.float64_t, ndim=2] buffer_out = np.zeros((sample, 3))
+    cdef np.float64_t[:, :] buffer = buffer_out
+
+    for i in range(sample):
+        a = np.random.rand()
+        b = np.random.rand()
+        c = np.random.rand()
+
+        while a + b + c > 1:
+            a = np.random.rand()
+            b = np.random.rand()
+            c = np.random.rand()
+
+        buffer[i, 0] = a
+        buffer[i, 1] = b
+        buffer[i, 2] = c
+
+    return buffer_out
+
+
+@cython.boundscheck(False)
+cdef int convex_hull(np.float64_t[4][2] ABCD, int[4] hull) nogil:
+    """Compute the convex hull of ABCD.
+    Parameters
+    ----------
+    ABCD : (4, 2) array_like
+        The four points
+    Returns
+    -------
+    hull : 2D array_like
+        The convex hull of ABCD
+    """
+    cdef int i, ilo
+    ilo = 0
+    for i in range(1, 4):
+        if ABCD[i][0] < ABCD[ilo][0]:
+            ilo = i
+    hull[0] = ilo
+    cdef np.uint8_t[4] mask
+    mask[0] = 1
+    mask[1] = 1
+    mask[2] = 1
+    mask[3] = 1
+    mask[hull[0]] = 0
+
+    # Iterate until we get back to the first point
+    cdef np.float64_t prev_angle = M_PI_2
+    cdef np.float64_t min_diff, diff, angle, best_angle
+    cdef int best_j
+
+    for i in range(1, 4):
+        min_diff = twopi
+        # Find point with smallest angle compared to previous segment
+        for j in range(4):
+            if j == hull[i-1]:  # Skip previous point
+                continue
+
+            angle = atan2(ABCD[j][1] - ABCD[hull[i-1]][1], ABCD[j][0] - ABCD[hull[i-1]][0])
+            diff = (prev_angle - angle) % twopi
+            if diff < min_diff:
+                min_diff = diff
+                best_j = j
+                best_angle = angle
+
+        if best_j == hull[0]:
+            break
+
+        hull[i] = best_j
+        mask[best_j] = 0
+        prev_angle = best_angle
+
+    for i in range(4):
+        if mask[i]:
+            hull[3] = i
+            return 3
+    return 4
+
+@cython.boundscheck(False)
+cdef int integrate_quad(
+    np.float64_t[2] A,
+    np.float64_t[2] B,
+    np.float64_t[2] C,
+    np.float64_t[2] P,
+    np.float64_t value,
+    np.float64_t[:, ::1] buffer,
+    np.uint8_t[:, ::1] buffer_mask,
+    int xs,
+    int ys,
+) nogil:
+    # Find bounding box of ABC
+    cdef int min_x = <int> (floor(xs * min(A[0], B[0], C[0])))
+    cdef int max_x = <int> (ceil(xs * max(A[0], B[0], C[0])))
+    cdef int min_y = <int> (floor(ys * min(A[1], B[1], C[1])))
+    cdef int max_y = <int> (ceil(ys * max(A[1], B[1], C[1])))
+
+    cdef np.float64_t[2] AB
+    cdef np.float64_t[2] CA
+    cdef np.float64_t[2] BC
+    AB[0] = B[0] - A[0]
+    AB[1] = B[1] - A[1]
+    CA[0] = A[0] - C[0]
+    CA[1] = A[1] - C[1]
+    BC[0] = C[0] - B[0]
+    BC[1] = C[1] - B[1]
+
+    cdef np.float64_t[2] PA
+    cdef np.float64_t[2] PB
+    cdef np.float64_t[2] PC
+    PA[0] = A[0] - P[0]
+    PA[1] = A[1] - P[1]
+    PB[0] = B[0] - P[0]
+    PB[1] = B[1] - P[1]
+    PC[0] = C[0] - P[0]
+    PC[1] = C[1] - P[1]
+
+    cdef np.float64_t[2] X
+    X[0] = 0
+    X[1] = 0
+
+    cdef int i, j
+    cdef np.float64_t[2] PX
+
+    cdef np.float64_t aa, bb, gg, alpha, beta, gamma, v, Vtot
+    Vtot = 0
+    # First pass: compute integral within ABC
+    for i in range(min_x, max_x):
+        X[0] = (i + 0.5) / xs
+        for j in range(min_y, max_y):
+            X[1] = (j + 0.5) / ys
+            PX[0] = X[0] - P[0]
+            PX[1] = X[1] - P[1]
+
+            # Compute alpha = (PX x BC) / (PB x BC)
+            if (PB[1] * BC[0] - PB[0] * BC[1]) == 0:
+                continue  # raise ValueError("PB x BC is zero")
+            alpha = 1 - (PX[1] * BC[0] - PX[0] * BC[1]) / (PB[1] * BC[0] - PB[0] * BC[1])
+            # Compute beta = (PX x CA) / (PC x CA)
+            if (PC[1] * CA[0] - PC[0] * CA[1]) == 0:
+                continue  # raise ValueError("PC x CA is zero")
+            beta = 1 - (PX[1] * CA[0] - PX[0] * CA[1]) / (PC[1] * CA[0] - PC[0] * CA[1])
+            # Compute gamma = (PX x AB) / (PA x AB)
+            if (PA[1] * AB[0] - PA[0] * AB[1]) == 0:
+                continue  # raise ValueError("PA x AB is zero")
+            gamma = 1 - (PX[1] * AB[0] - PX[0] * AB[1]) / (PA[1] * AB[0] - PA[0] * AB[1])
+
+            # Check whether X is within ABC
+            aa = 0 <= alpha
+            bb = 0 <= beta
+            gg = 0 <= gamma
+
+            if not (aa and bb and gg):
+                continue
+
+            v = 2
+            if 0 <= alpha < 1:
+                v = alpha
+            if 0 <= beta < 1:
+                v = min(v, beta)
+            if 0 <= gamma < 1:
+                v = min(v, gamma)
+
+            Vtot += v
+
+    # Special case: Vtot == 0
+    if Vtot == 0:
+        _add_point_to_greyscale_image(
+            buffer, buffer_mask, P[0], P[1], value, xs, ys
+        )
+        return 0
+
+    # Second pass: deposit
+    for i in range(min_x, max_x):
+        X[0] = (i + 0.5) / xs
+        for j in range(min_y, max_y):
+            X[1] = (j + 0.5) / ys
+            PX[0] = X[0] - P[0]
+            PX[1] = X[1] - P[1]
+
+            # Compute alpha = (PX x BC) / (PB x BC)
+            alpha = 1 - (PX[1] * BC[0] - PX[0] * BC[1]) / (PB[1] * BC[0] - PB[0] * BC[1])
+            # Compute beta = (PX x CA) / (PC x CA)
+            beta = 1 - (PX[1] * CA[0] - PX[0] * CA[1]) / (PC[1] * CA[0] - PC[0] * CA[1])
+            # Compute gamma = (PX x AB) / (PA x AB)
+            gamma = 1 - (PX[1] * AB[0] - PX[0] * AB[1]) / (PA[1] * AB[0] - PA[0] * AB[1])
+
+            # Check whether X is within ABC
+            aa = 0 <= alpha
+            bb = 0 <= beta
+            gg = 0 <= gamma
+
+            if not (aa and bb and gg):
+                continue
+
+            v = 2
+            if 0 <= alpha < 1:
+                v = alpha
+            if 0 <= beta < 1:
+                v = min(v, beta)
+            if 0 <= gamma < 1:
+                v = min(v, gamma)
+
+            buffer[j % ys, i % xs] += value * v / Vtot
+            buffer_mask[j % ys, i % xs] = 1
+
+@cython.boundscheck(False)
+cdef int integrate_tetrahedron_proj(
+    np.float64_t[3] origin,
+    np.float64_t[3] a,
+    np.float64_t[3] b,
+    np.float64_t[3] c,
+    np.float64_t value,
+    np.float64_t[:, ::1] buffer,
+    np.uint8_t[:, ::1] buffer_mask,
+    int xs,
+    int ys
+) nogil:
+    cdef np.float64_t[4][2] ABCD_xy
+    cdef np.float64_t[2] A_xy, B_xy, C_xy, D_xy, P_xy, P1_xy, P2_xy
+    cdef np.float64_t det, S1, S2
+
+    ABCD_xy[0][0] = origin[0]
+    ABCD_xy[0][1] = origin[1]
+    ABCD_xy[1][0] = origin[0] + a[0]
+    ABCD_xy[1][1] = origin[1] + a[1]
+    ABCD_xy[2][0] = origin[0] + b[0]
+    ABCD_xy[2][1] = origin[1] + b[1]
+    ABCD_xy[3][0] = origin[0] + c[0]
+    ABCD_xy[3][1] = origin[1] + c[1]
+    P_xy[0] = 0
+    P_xy[1] = 0
+
+    cdef int[4] iABC
+
+    cdef int Nhull = convex_hull(ABCD_xy, iABC)
+
+    A_xy[0] = ABCD_xy[iABC[0]][0]
+    A_xy[1] = ABCD_xy[iABC[0]][1]
+    B_xy[0] = ABCD_xy[iABC[1]][0]
+    B_xy[1] = ABCD_xy[iABC[1]][1]
+    C_xy[0] = ABCD_xy[iABC[2]][0]
+    C_xy[1] = ABCD_xy[iABC[2]][1]
+    D_xy[0] = ABCD_xy[iABC[3]][0]
+    D_xy[1] = ABCD_xy[iABC[3]][1]
+    if Nhull == 4:
+        # Find the intersection of AC with BD
+        det = ((A_xy[0] - C_xy[0]) * (B_xy[1] - D_xy[1]) - (A_xy[1] - C_xy[1]) * (B_xy[0] - D_xy[0]))
+        P_xy[0] = ((A_xy[0] * C_xy[1] - A_xy[1] * C_xy[0]) * (B_xy[0] - D_xy[0]) - (A_xy[0] - C_xy[0]) * (B_xy[0] * D_xy[1] - B_xy[1] * D_xy[0])) / det
+        P_xy[1] = ((A_xy[0] * C_xy[1] - A_xy[1] * C_xy[0]) * (B_xy[1] - D_xy[1]) - (A_xy[1] - C_xy[1]) * (B_xy[0] * D_xy[1] - B_xy[1] * D_xy[0])) / det
+
+        # Compute (double the) surface of each triangle
+        # Surface of triangle ABC
+        S1 = (C_xy[1] - A_xy[1]) * (B_xy[0] - A_xy[0]) - (C_xy[0] - A_xy[0]) * (B_xy[1] - A_xy[1])
+        # Surface of triangle ACD
+        S2 = (D_xy[1] - A_xy[1]) * (C_xy[0] - A_xy[0]) - (D_xy[0] - A_xy[0]) * (C_xy[1] - A_xy[1])
+
+        # Move slightly towards B to prevent rounding errors
+        P1_xy[0] = P_xy[0] + 1e-10 * (B_xy[0] - P_xy[0])
+        P1_xy[1] = P_xy[1] + 1e-10 * (B_xy[1] - P_xy[1])
+        integrate_quad(A_xy, B_xy, C_xy, P1_xy, value * S1 / (S1 + S2), buffer, buffer_mask, xs, ys)
+
+        # Move slightly towards D to prevent rounding errors
+        P2_xy[0] = P_xy[0] + 1e-10 * (D_xy[0] - P_xy[0])
+        P2_xy[1] = P_xy[1] + 1e-10 * (D_xy[1] - P_xy[1])
+        integrate_quad(A_xy, C_xy, D_xy, P2_xy, value * S2 / (S1 + S2), buffer, buffer_mask, xs, ys)
+    else:
+        integrate_quad(A_xy, B_xy, C_xy, D_xy, value, buffer, buffer_mask, xs, ys)
+
+# @cython.boundscheck(False)
+def add_points_to_greyscale_image_with_lagrangian_tesselation(
+    np.ndarray[np.float64_t, ndim=2] buffer,
+    np.ndarray[np.uint8_t,   ndim=2] buffer_mask,
+    np.ndarray[np.float64_t, ndim=4] p3d,
+    np.ndarray[np.float64_t, ndim=3] pv,
+):
+    cdef int xs = buffer.shape[0]
+    cdef int ys = buffer.shape[1]
+    cdef int Ngrid = p3d.shape[0]
+
+    cdef np.float64_t[:, ::1] buffer_view = buffer
+    cdef np.uint8_t[:, ::1] buffer_mask_view = buffer_mask
+    cdef np.float64_t[:, :, :, ::1] p3d_view = p3d
+
+    cdef np.ndarray[np.uint8_t, ndim=2] vert = np.array(
+        (
+            (0, 0, 0),
+            (1, 0, 0),
+            (1, 1, 0),
+            (0, 1, 0),
+            (0, 0, 1),
+            (1, 0, 1),
+            (1, 1, 1),
+            (0, 1, 1),
+        ),
+        dtype=np.uint8
+    )
+    cdef np.ndarray[np.uint8_t, ndim=2] conn = np.array(
+        (
+            (4, 0, 7, 1),
+            (1, 0, 3, 7),
+            (5, 1, 4, 7),
+            (2, 3, 7, 1),
+            (1, 5, 6, 7),
+            (2, 6, 1, 7),
+        ),
+        dtype=np.uint8
+    )
+
+    cdef int i, m
+    cdef np.uint8_t[:, ::1] off
+
+    for m in range(len(conn)):
+        # ro = sample_tetrahedron(split)
+        off = vert[conn[m]]
+
+        for i in prange(Ngrid, nogil=True):
+            # print(f"{m=}/{len(conn)} {i=}/{Ngrid}")
+            _lagrangian_tesselation_helper(
+                Ngrid,
+                i,
+                p3d_view,
+                off,
+                buffer_view,
+                buffer_mask_view,
+                xs,
+                ys
+            )
+
+
+@cython.boundscheck(False)
+cdef int _lagrangian_tesselation_helper(
+    const int Ngrid,
+    const int i,
+    const np.float64_t[:, :, :, ::1] p3d_view,
+    const np.uint8_t[:, ::1] off,
+    np.float64_t[:, ::1] buffer_view,
+    np.uint8_t[:, ::1] buffer_mask_view,
+    const int xs,
+    const int ys,
+) nogil:
+    cdef int j, k, idim
+    cdef np.float64_t[3] orig
+    cdef np.float64_t[3] a
+    cdef np.float64_t[3] b
+    cdef np.float64_t[3] c
+
+    for j in range(Ngrid):
+        for k in range(Ngrid):
+            for idim in range(3):   # stop at 2, because z is not used
+                orig[idim] = p3d_view[(i + off[3][0]) % Ngrid, (j + off[3][1]) % Ngrid, (k + off[3][2]) % Ngrid, idim]
+                a[idim] = _wrap_dist(p3d_view[(i + off[0][0]) % Ngrid, (j + off[0][1]) % Ngrid, (k + off[0][2]) % Ngrid, idim] - orig[idim])
+                b[idim] = _wrap_dist(p3d_view[(i + off[1][0]) % Ngrid, (j + off[1][1]) % Ngrid, (k + off[1][2]) % Ngrid, idim] - orig[idim])
+                c[idim] = _wrap_dist(p3d_view[(i + off[2][0]) % Ngrid, (j + off[2][1]) % Ngrid, (k + off[2][2]) % Ngrid, idim] - orig[idim])
+
+            integrate_tetrahedron_proj(
+                orig,
+                a,
+                b,
+                c,
+                1,
+                buffer_view,
+                buffer_mask_view,
+                xs,
+                ys
+            )
diff --git a/yt/visualization/fixed_resolution.py b/yt/visualization/fixed_resolution.py
index d2da6a29e1..9adf82909d 100644
--- a/yt/visualization/fixed_resolution.py
+++ b/yt/visualization/fixed_resolution.py
@@ -11,9 +11,11 @@
 from yt.frontends.ytdata.utilities import save_as_dataset
 from yt.funcs import get_output_filename, iter_fields, mylog
 from yt.loaders import load_uniform_grid
+from yt.utilities.exceptions import YTException
 from yt.utilities.lib.api import (  # type: ignore
     CICDeposit_2,
     add_points_to_greyscale_image,
+    add_points_to_greyscale_image_with_lagrangian_tesselation,
 )
 from yt.utilities.lib.pixelization_routines import (
     pixelize_cylinder,
@@ -719,9 +721,11 @@ def __init__(
         if axis is not None:
             self.xax = self.ds.coordinates.x_axis[axis]
             self.yax = self.ds.coordinates.y_axis[axis]
+            self.zax = list({0, 1, 2} - {self.xax, self.yax})[0]
             axis_name = self.ds.coordinates.axis_name
             self.x_field = f"particle_position_{axis_name[self.xax]}"
             self.y_field = f"particle_position_{axis_name[self.yax]}"
+            self.z_field = f"particle_position_{axis_name[self.zax]}"
 
     @override
     def _generate_image_and_mask(self, item) -> None:
@@ -810,6 +814,43 @@ def _generate_image_and_mask(self, item) -> None:
                 x_bin_edges,
                 y_bin_edges,
             )
+        elif deposition == "lagrangian_tesselation":
+            z_data = self.data_source.dd[ftype, self.z_field]
+            dz = z_data.in_units("code_length").d
+            # TODO: handle periodicity
+            pz = dz / (bounds[1] - bounds[0])
+            order = np.argsort(self.data_source.dd[ftype, "particle_index"])
+
+            # Check all particles have the same mass
+            if not np.allclose(
+                self.data_source.dd[ftype, "particle_mass"].v,
+                self.data_source.dd[ftype, "particle_mass"][0].v,
+            ):
+                raise YTException(
+                    "Particles must have the same mass for Lagrangian "
+                    "tesselation to make sense."
+                )
+
+            # Check particles are compatible with being initially on a grid
+            Ngrid = 1
+            while Ngrid**3 < len(order):
+                Ngrid += 1
+
+            if Ngrid**3 != len(order):
+                raise YTException(
+                    "Particles must initially reside on a grid for Lagrangian "
+                    "tesselation to make sense."
+                )
+            p3d = np.stack([px[order], py[order], pz[order]], axis=1).reshape(
+                Ngrid, Ngrid, Ngrid, 3
+            )
+
+            add_points_to_greyscale_image_with_lagrangian_tesselation(
+                buff,
+                buff_mask,
+                p3d.copy(),
+                splat_vals[None, None, :],
+            )
         else:
             raise ValueError(f"Received unknown deposition method '{deposition}'")