diff --git a/src/quantem/core/config.py b/src/quantem/core/config.py
index c2eda10b..7213f3b5 100644
--- a/src/quantem/core/config.py
+++ b/src/quantem/core/config.py
@@ -34,6 +34,15 @@
     if "cuda" in str(e):
         NUM_DEVICES = 0
     _defaults["has_cupy"] = False
+try:
+    import quantem.cuda  # type: ignore  # noqa: F401
+
+    _defaults["has_quantem_cuda"] = True
+except ModuleNotFoundError:
+    _defaults["has_quantem_cuda"] = False
+except Exception:
+    # installed but unloadable (e.g. libcudart missing at runtime)
+    _defaults["has_quantem_cuda"] = False
 
 
 defaults: list[Mapping] = [_defaults]
diff --git a/src/quantem/core/ml/models/kplanes.py b/src/quantem/core/ml/models/kplanes.py
index cdf55261..96ab9dd3 100644
--- a/src/quantem/core/ml/models/kplanes.py
+++ b/src/quantem/core/ml/models/kplanes.py
@@ -10,6 +10,8 @@
 import torch.nn.functional as F
 from torch import nn
 
+from quantem.core import config
+
 from .model_base import PPLR, TensorDecompositionModel
 from .so3params import SO3ParamQuat, SO3ParamR9SVD
 
@@ -169,6 +171,7 @@ class KPlanes(PPLR, TensorDecompositionModel):
     """
     K-Planes model adapted from Fridovich-Keil et al., https://arxiv.org/abs/2301.10241
     """
+
     def __init__(
         self,
         # Grid parameters
@@ -312,10 +315,30 @@ def interpolate_ms_features_tilted(
     """
     Fully-vectorized multi-scale, multi-rotation K-Planes feature interpolation.
     Returns features of shape (B, C * T * num_scales).
+
+    When the optional ``quantem-cuda`` package is installed, the tensors live
+    on a CUDA device, and the ``use_cuda_kernels`` config option is true
+    (default), each scale dispatches to the fused CUDA kernel (rotate +
+    grid_sample + Hadamard product in one launch, analytic backward) —
+    identical semantics to the torch path below.
     """
     T = rotation_matrices.shape[0]
     B = pts.shape[0]
 
+    if (
+        pts.is_cuda
+        and pts.dtype == torch.float32
+        and rotation_matrices.dtype == torch.float32
+        and all(g.dtype == torch.float32 for g in ms_grids)
+        and config.get("has_quantem_cuda")
+        and config.get("use_cuda_kernels", default=True)
+    ):
+        from quantem.cuda.core.ml import kplanes_tilted_fuse
+
+        return torch.cat(
+            [kplanes_tilted_fuse(pts, rotation_matrices, g) for g in ms_grids], dim=-1
+        )
+
     # (T, B, 3)  — rotate all points by all rotations at once
     rotated = torch.einsum("tij,bj->tbi", rotation_matrices, pts)
 
diff --git a/src/quantem/core/quantem.yaml b/src/quantem/core/quantem.yaml
index 80f34587..0e4913a4 100644
--- a/src/quantem/core/quantem.yaml
+++ b/src/quantem/core/quantem.yaml
@@ -14,6 +14,10 @@ dtype_complex: complex64
 # global verbosity, not currently used
 
 verbose: 1
+
+# Dispatch to the fused quantem-cuda kernels when installed (pip install quantem[cuda])
+# and the tensors are on a CUDA device. Set false to force the pure-torch paths.
+use_cuda_kernels: true
 cupy:
   # The size of the fft cache in MB used by cupy
   # https://docs.cupy.dev/en/stable/user_guide/fft.html#fft-plan-cache
diff --git a/src/quantem/diffractive_imaging/object_models.py b/src/quantem/diffractive_imaging/object_models.py
index d311d8dd..17c4e72c 100644
--- a/src/quantem/diffractive_imaging/object_models.py
+++ b/src/quantem/diffractive_imaging/object_models.py
@@ -396,6 +396,30 @@ def get_tv_loss(
         return loss
 
     def _calc_tv_loss(self, array: torch.Tensor, weight: tuple[float, float]) -> torch.Tensor:
+        # Identical math via the fused quantem-cuda L1 kernel when available.
+        # Weighted size-1 axes fall through to torch so degenerate inputs
+        # behave exactly as before.
+        if (
+            array.ndim == 3
+            and array.is_cuda
+            and array.dtype == torch.float32
+            and config.get("has_quantem_cuda")
+            and config.get("use_cuda_kernels", default=True)
+            and not (weight[0] > 0 and array.shape[0] == 1)
+            and not (weight[1] > 0 and (array.shape[1] == 1 or array.shape[2] == 1))
+        ):
+            from quantem.cuda.core import tv_loss_l1_3d
+
+            s, h, w_ = array.shape
+            wts = array.new_tensor([weight[0], weight[1], weight[1]])
+            active = int((wts > 0).sum())
+            if active == 0:
+                return self._get_zero_loss_tensor()
+            counts = array.new_tensor(
+                [(s - 1) * h * w_, s * (h - 1) * w_, s * h * (w_ - 1)]
+            ).clamp(min=1)
+            return (wts * tv_loss_l1_3d(array) / counts).sum() / active
+
         loss = self._get_zero_loss_tensor()
         calc_dim = 0
         for dim in range(array.ndim):
diff --git a/src/quantem/tomography/object_models.py b/src/quantem/tomography/object_models.py
index 817e7ed5..41e9793e 100644
--- a/src/quantem/tomography/object_models.py
+++ b/src/quantem/tomography/object_models.py
@@ -18,6 +18,7 @@
 from quantem.core.utils.rng import RNGMixin
 from quantem.tomography.dataset_models import TomographyINRPretrainDataset
 from quantem.tomography.tomography_context import ReconstructionContext
+from quantem.tomography.utils import tv_loss_vol_sq
 
 
 class ObjConstraintParams:
@@ -439,11 +440,9 @@ def get_tv_loss(self, ctx: ReconstructionContext) -> torch.Tensor:
         # TV over the three trailing spatial dims, leaving any leading channel/batch axes
         # intact. Works for a 3-D volume, obj_view's [1, D, H, W], and a multimodal
         # [C, D, H, W] (channels = elemental compositions), matching the INR / tensor-decomp
-        # convention where the object carries a leading channel dimension.
-        tv_d = torch.pow(ctx.obj[..., 1:, :, :] - ctx.obj[..., :-1, :, :], 2).sum()
-        tv_h = torch.pow(ctx.obj[..., :, 1:, :] - ctx.obj[..., :, :-1, :], 2).sum()
-        tv_w = torch.pow(ctx.obj[..., :, :, 1:] - ctx.obj[..., :, :, :-1], 2).sum()
-        tv_loss = tv_d + tv_h + tv_w
+        # convention where the object carries a leading channel dimension. tv_loss_vol_sq
+        # dispatches to the fused quantem-cuda kernel when available.
+        tv_loss = tv_loss_vol_sq(ctx.obj)
 
         return tv_loss * self.constraints.tv_vol / ctx.obj.numel()
 
@@ -573,6 +572,10 @@ def apply_hard_constraints(self, pred: torch.Tensor) -> torch.Tensor:
 
         return pred
 
+    def sample_tv_tap_coords(self, coords: torch.Tensor) -> Optional[torch.Tensor]:
+        """Hook for the training loop: returns None (INR TV uses autograd, no tap merging)."""
+        return None
+
     # --- Define get_tv_loss ---
 
     def get_tv_loss(self, ctx: ReconstructionContext) -> torch.Tensor:
@@ -674,6 +677,35 @@ def forward(self, coords: Optional[torch.Tensor] = None) -> torch.Tensor:
 
         return all_densities
 
+    def forward_with_tv_taps(
+        self, coords: torch.Tensor, tap_coords: torch.Tensor
+    ) -> tuple[torch.Tensor, torch.Tensor]:
+        """
+        Single model call covering the main batch and the volume-TV tap points.
+
+        The out-of-bounds mask and hard constraints apply to the main chunk
+        only; tap densities are returned raw (border-clamped), matching the
+        fallback path in ``get_volume_tv_loss``.
+        """
+        merged = self.model(torch.cat([coords, tap_coords], dim=0))
+        if isinstance(merged, tuple):
+            merged = merged[0]
+        main, taps = merged[: coords.shape[0]], merged[coords.shape[0] :]
+
+        if main.dim() > 1:
+            main = main.squeeze(-1)
+        valid_mask = (
+            (coords[:, 0] >= -1) & (coords[:, 0] <= 1) & (coords[:, 1] >= -1) & (coords[:, 1] <= 1)
+        ).float()
+        if main.dim() > 1:
+            valid_mask = valid_mask.unsqueeze(-1)
+        main = main * valid_mask
+        main = self.apply_hard_constraints(main)
+
+        if taps.dim() == 1:
+            taps = taps.unsqueeze(-1)
+        return main, taps
+
     # Pretrain Loop
 
     def pretrain(
@@ -911,13 +943,37 @@ def from_model(
         obj_model.to(device)
         return obj_model
 
+    def sample_tv_tap_coords(self, coords: torch.Tensor) -> Optional[torch.Tensor]:
+        """
+        Sample the finite-difference tap coordinates for the volume TV loss.
+
+        Returns a (4*n, 3) tensor [base; base+h*ex; base+h*ey; base+h*ez] for n
+        sampled base points, or None when tv_vol == 0.  The training loop
+        concatenates this to all_coords so the TV taps are evaluated in the
+        same model call as the main forward pass.
+        """
+        if self.constraints.tv_vol == 0:
+            return None
+        model = _unwrap(self.model)
+        h = 2.0 / min(model.resolution)
+        num_tv_samples = min(10_000, coords.shape[0])
+        tv_indices = torch.randperm(coords.shape[0], device=coords.device)[:num_tv_samples]
+        tv_coords = coords[tv_indices]  # (n, 3)
+        ex = torch.zeros(3, device=tv_coords.device)
+        ex[0] = h
+        ey = torch.zeros(3, device=tv_coords.device)
+        ey[1] = h
+        ez = torch.zeros(3, device=tv_coords.device)
+        ez[2] = h
+        return torch.cat([tv_coords, tv_coords + ex, tv_coords + ey, tv_coords + ez], dim=0)
+
     # --- Constraints ---
 
     def apply_soft_constraints(self, ctx: ReconstructionContext) -> torch.Tensor:
         soft_loss = torch.tensor(
             0.0, device=ctx.pred.device if ctx.pred is not None else self.device
         )
-        if self.constraints.tv_vol > 0:
+        if self.constraints.tv_plane > 0 or self.constraints.tv_vol > 0:
             assert ctx.coords is not None, "Coordinates must be provided for TV loss"
             assert ctx.pred is not None, "Prediction must be provided for TV loss"
             soft_loss += self.get_tv_loss(ctx)
@@ -943,8 +999,12 @@ def get_tv_loss(self, ctx: ReconstructionContext) -> torch.Tensor:
         assert ctx.coords is not None, "Coordinates must be provided for TV loss"
         assert ctx.pred is not None, "Prediction must be provided for TV loss"
         tv_loss = torch.tensor(0.0, device=ctx.pred.device)
-        tv_loss += self._get_plane_tv_loss()
-        tv_loss += self.get_volume_tv_loss(ctx.coords)
+        if self.constraints.tv_plane > 0:
+            tv_loss += self._get_plane_tv_loss()
+        if self.constraints.tv_vol > 0:
+            tv_loss += self.get_volume_tv_loss(
+                ctx.coords, precomputed_tap_densities=ctx.tv_tap_densities
+            )
         return tv_loss
 
     def _get_plane_tv_loss(self) -> torch.Tensor:
@@ -972,37 +1032,52 @@ def _get_plane_tv_loss(self) -> torch.Tensor:
 
         return self.constraints.tv_plane * torch.stack(per_level).sum()
 
-    def get_volume_tv_loss(self, coords: torch.Tensor) -> torch.Tensor:
+    def get_volume_tv_loss(
+        self,
+        coords: torch.Tensor,
+        precomputed_tap_densities: Optional[torch.Tensor] = None,
+    ) -> torch.Tensor:
         """
         Isotropic volume TV via finite differences. Same form as the autograd
         version (L1 of gradient L2-norm) but avoids double-backward, so it
         works for KPlanesTILTED, CPTilted, and anything else.
-        """
-        num_tv_samples = min(10_000, coords.shape[0])
-        tv_indices = torch.randperm(coords.shape[0], device=coords.device)[:num_tv_samples]
-        tv_coords = coords[tv_indices]  # (N, 3)
 
+        When *precomputed_tap_densities* is provided (a (4N, C) tensor from the
+        merged single-pass forward in the training loop), the model call is
+        skipped entirely and the supplied values are used directly.  When absent
+        the existing batched 4N-point fallback path runs unchanged.
+        """
         model = _unwrap(self.model)
         h = 2.0 / min(model.resolution)
 
-        pred = model(tv_coords)
-        if isinstance(pred, tuple):
-            pred = pred[0]
-        if pred.dim() == 1:
-            pred = pred.unsqueeze(-1)  # (N, 1)
-
-        grads = []
-        for axis in range(3):
-            offset = torch.zeros(3, device=tv_coords.device)
-            offset[axis] = h
-            shifted_pred = self.model(tv_coords + offset)
-            if isinstance(shifted_pred, tuple):
-                shifted_pred = shifted_pred[0]
-            if shifted_pred.dim() == 1:
-                shifted_pred = shifted_pred.unsqueeze(-1)
-            grads.append((shifted_pred - pred) / h)  # (N, 1)
-
-        grad_stack = torch.stack(grads, dim=-1)  # (N, C, 3)
+        if precomputed_tap_densities is not None:
+            batched_pred = precomputed_tap_densities
+        else:
+            num_tv_samples = min(10_000, coords.shape[0])
+            tv_indices = torch.randperm(coords.shape[0], device=coords.device)[:num_tv_samples]
+            tv_coords = coords[tv_indices]  # (N, 3)
+
+            ex = torch.zeros(3, device=tv_coords.device)
+            ex[0] = h
+            ey = torch.zeros(3, device=tv_coords.device)
+            ey[1] = h
+            ez = torch.zeros(3, device=tv_coords.device)
+            ez[2] = h
+
+            batched = torch.cat(
+                [tv_coords, tv_coords + ex, tv_coords + ey, tv_coords + ez], dim=0
+            )  # (4N, 3)
+            batched_pred = model(batched)
+            if isinstance(batched_pred, tuple):
+                batched_pred = batched_pred[0]
+
+        if batched_pred.dim() == 1:
+            batched_pred = batched_pred.unsqueeze(-1)  # (4N, C)
+
+        pred, px, py, pz = batched_pred.chunk(4, dim=0)  # each (N, C)
+        grad_stack = torch.stack(
+            [(px - pred) / h, (py - pred) / h, (pz - pred) / h], dim=-1
+        )  # (N, C, 3)
         grad_norm = torch.norm(grad_stack, dim=-1)  # (N, C)
 
         return self.constraints.tv_vol * grad_norm.mean()
diff --git a/src/quantem/tomography/tomography.py b/src/quantem/tomography/tomography.py
index f86095ad..a3789f63 100644
--- a/src/quantem/tomography/tomography.py
+++ b/src/quantem/tomography/tomography.py
@@ -210,7 +210,14 @@ def reconstruct(
                 ):
                     all_coords = self.dset.get_coords(batch, N, curr_num_samples_per_ray)
 
-                    all_densities = self.obj_model.forward(all_coords)
+                    tap_coords = self.obj_model.sample_tv_tap_coords(all_coords)
+                    if tap_coords is not None:
+                        all_densities, tv_tap_raw = self.obj_model.forward_with_tv_taps(
+                            all_coords, tap_coords
+                        )
+                    else:
+                        all_densities = self.obj_model.forward(all_coords)
+                        tv_tap_raw = None
 
                     integrated_densities = self.dset.integrate_rays(
                         all_densities,
@@ -225,6 +232,7 @@ def reconstruct(
                         coords=all_coords,
                         pred=pred,
                         all_densities=all_densities,
+                        tv_tap_densities=tv_tap_raw,
                     )
                 )
 
diff --git a/src/quantem/tomography/tomography_context.py b/src/quantem/tomography/tomography_context.py
index d322287c..56a9b262 100644
--- a/src/quantem/tomography/tomography_context.py
+++ b/src/quantem/tomography/tomography_context.py
@@ -1,9 +1,10 @@
 from dataclasses import dataclass
 from typing import Optional
-from quantem.core.ml.constraints import BaseContext
 
 import torch
 
+from quantem.core.ml.constraints import BaseContext
+
 
 @dataclass
 class ReconstructionContext(BaseContext):
@@ -28,3 +29,4 @@ class ReconstructionContext(BaseContext):
     pred: Optional[torch.Tensor] = None
     all_densities: Optional[torch.Tensor] = None
     obj: Optional[torch.Tensor] = None
+    tv_tap_densities: Optional[torch.Tensor] = None
diff --git a/src/quantem/tomography/utils.py b/src/quantem/tomography/utils.py
index 08093925..6c8d3491 100644
--- a/src/quantem/tomography/utils.py
+++ b/src/quantem/tomography/utils.py
@@ -1,6 +1,8 @@
 import torch
 import torch.nn.functional as F
 
+from quantem.core import config
+
 # --- Projection Operator Utils ---
 
 
@@ -69,6 +71,42 @@ def transform_slice(mag_slice):
     return rotated_mags.permute(1, 2, 3, 0)
 
 
+def tv_loss_vol_sq(obj: torch.Tensor) -> torch.Tensor:
+    """Squared-anisotropic volume TV: sum of squared forward differences.
+
+    Computes ``Σ (Δd)² + Σ (Δh)² + Σ (Δw)²`` over the three trailing
+    spatial dims, leaving any leading channel/batch axes intact (they are
+    included in the sum). This is the unnormalized ``tv_vol`` regularizer;
+    callers apply their own ``weight / numel`` scaling.
+
+    When the optional ``quantem-cuda`` package is installed
+    (``pip install quantem[cuda]``), the tensor is on a CUDA device, and the
+    ``use_cuda_kernels`` config option is true (default), this dispatches to
+    the fused CUDA forward/backward kernel — identical math, one kernel
+    launch instead of several large intermediates.
+
+    Args:
+        obj: Tensor of shape ``[..., D, H, W]`` (ndim >= 3).
+
+    Returns:
+        0-dim tensor on the same device as ``obj``; differentiable.
+    """
+    if (
+        obj.is_cuda
+        and obj.dtype == torch.float32
+        and config.get("has_quantem_cuda")
+        and config.get("use_cuda_kernels", default=True)
+    ):
+        from quantem.cuda.core import tv_loss_sq_3d
+
+        return tv_loss_sq_3d(obj)
+
+    tv_d = torch.pow(obj[..., 1:, :, :] - obj[..., :-1, :, :], 2).sum()
+    tv_h = torch.pow(obj[..., :, 1:, :] - obj[..., :, :-1, :], 2).sum()
+    tv_w = torch.pow(obj[..., :, :, 1:] - obj[..., :, :, :-1], 2).sum()
+    return tv_d + tv_h + tv_w
+
+
 def tv_loss_1d(x: torch.Tensor, reduction: str = "mean") -> torch.Tensor:
     """
     1D Total Variation Loss.
diff --git a/tests/diffractive_imaging/test_tv_cuda_dispatch.py b/tests/diffractive_imaging/test_tv_cuda_dispatch.py
new file mode 100644
index 00000000..e320ceb3
--- /dev/null
+++ b/tests/diffractive_imaging/test_tv_cuda_dispatch.py
@@ -0,0 +1,103 @@
+"""Tests for ``ObjectConstraints._calc_tv_loss``'s optional quantem-cuda dispatch.
+
+The constraint must produce the same loss and gradients on every path: pure
+torch on CPU, pure torch on GPU (kill-switch or quantem-cuda absent), and
+the fused L1 kernel when quantem-cuda is installed. The dispatch itself is
+asserted by monkeypatching the kernel entry point, so these tests are
+meaningful both with and without quantem-cuda in the environment.
+"""
+
+import pytest
+import torch
+
+from quantem.core import config
+from quantem.diffractive_imaging.object_models import ObjectConstraints
+
+requires_gpu = pytest.mark.skipif(not torch.cuda.is_available(), reason="requires a CUDA device")
+requires_quantem_cuda = pytest.mark.skipif(
+    not config.get("has_quantem_cuda"), reason="requires quantem-cuda"
+)
+
+
+class _Host:
+    """Minimal stand-in providing what _calc_tv_loss needs from self."""
+
+    def __init__(self, device):
+        self.device = device
+
+    _get_zero_loss_tensor = ObjectConstraints._get_zero_loss_tensor
+    _calc_tv_loss = ObjectConstraints._calc_tv_loss
+
+
+def _phase(device, shape=(8, 24, 20), seed=0, requires_grad=False):
+    gen = torch.Generator().manual_seed(seed)
+    arr = torch.rand(shape, generator=gen, dtype=torch.float32)
+    return arr.to(device).requires_grad_(requires_grad)
+
+
+WEIGHTS = [(5.0, 0.1), (5.0, 0.0), (0.0, 0.3), (0.0, 0.0)]
+
+
+@requires_gpu
+@pytest.mark.parametrize("weight", WEIGHTS)
+def test_gpu_matches_cpu_reference(weight):
+    expected = _Host("cpu")._calc_tv_loss(_phase("cpu"), weight)
+    actual = _Host("cuda")._calc_tv_loss(_phase("cuda"), weight)
+    torch.testing.assert_close(actual.cpu(), expected, rtol=1e-5, atol=1e-7)
+
+
+@requires_gpu
+@pytest.mark.parametrize("weight", WEIGHTS[:3])
+def test_gpu_grads_match_cpu_reference(weight):
+    a_cpu = _phase("cpu", requires_grad=True)
+    _Host("cpu")._calc_tv_loss(a_cpu, weight).backward()
+    a_gpu = _phase("cuda", requires_grad=True)
+    _Host("cuda")._calc_tv_loss(a_gpu, weight).backward()
+    torch.testing.assert_close(a_gpu.grad.cpu(), a_cpu.grad, rtol=1e-4, atol=1e-7)
+
+
+@requires_gpu
+@pytest.mark.parametrize("shape", [(1, 24, 20), (24, 20)])
+def test_degenerate_and_2d_match_cpu(shape):
+    # num_slices == 1 (weight[0] zeroed upstream, as get_tv_loss does) and
+    # plain 2-D arrays keep exact parity whichever path runs.
+    weight = (0.0, 0.3)
+    expected = _Host("cpu")._calc_tv_loss(_phase("cpu", shape=shape), weight)
+    actual = _Host("cuda")._calc_tv_loss(_phase("cuda", shape=shape), weight)
+    torch.testing.assert_close(actual.cpu(), expected, rtol=1e-5, atol=1e-7)
+
+
+@requires_gpu
+@requires_quantem_cuda
+def test_dispatches_to_kernel(monkeypatch):
+    import quantem.cuda.core
+
+    calls = []
+    real = quantem.cuda.core.tv_loss_l1_3d
+
+    def spy(volume):
+        calls.append(volume.shape)
+        return real(volume)
+
+    monkeypatch.setattr(quantem.cuda.core, "tv_loss_l1_3d", spy)
+    _Host("cuda")._calc_tv_loss(_phase("cuda"), (5.0, 0.1))
+    assert len(calls) == 1
+
+
+@requires_gpu
+@requires_quantem_cuda
+def test_kill_switch_forces_torch_path(monkeypatch):
+    import quantem.cuda.core
+
+    def boom(volume):
+        raise AssertionError("kernel should not be called with use_cuda_kernels=False")
+
+    monkeypatch.setattr(quantem.cuda.core, "tv_loss_l1_3d", boom)
+    expected = _Host("cpu")._calc_tv_loss(_phase("cpu"), (5.0, 0.1))
+    # config.set lacks __exit__, so restore explicitly rather than via `with`.
+    config.set({"use_cuda_kernels": False})
+    try:
+        actual = _Host("cuda")._calc_tv_loss(_phase("cuda"), (5.0, 0.1))
+    finally:
+        config.set({"use_cuda_kernels": True})
+    torch.testing.assert_close(actual.cpu(), expected, rtol=1e-5, atol=1e-7)
diff --git a/tests/ml/test_kplanes_cuda_dispatch.py b/tests/ml/test_kplanes_cuda_dispatch.py
new file mode 100644
index 00000000..da61f6db
--- /dev/null
+++ b/tests/ml/test_kplanes_cuda_dispatch.py
@@ -0,0 +1,106 @@
+"""Tests for ``interpolate_ms_features_tilted``'s optional quantem-cuda dispatch.
+
+The function must produce the same features and gradients on every path:
+pure torch on CPU, pure torch on GPU (kill-switch or quantem-cuda absent),
+and the fused CUDA kernel when quantem-cuda is installed. The dispatch
+itself is asserted by monkeypatching the kernel entry point, so these tests
+are meaningful both with and without quantem-cuda in the environment.
+"""
+
+import pytest
+import torch
+from torch import nn
+
+from quantem.core import config
+from quantem.core.ml.models.kplanes import interpolate_ms_features_tilted
+
+requires_gpu = pytest.mark.skipif(not torch.cuda.is_available(), reason="requires a CUDA device")
+requires_quantem_cuda = pytest.mark.skipif(
+    not config.get("has_quantem_cuda"), reason="requires quantem-cuda"
+)
+
+
+def _inputs(device, requires_grad=False, seed=0, B=64, T=3, C=5, scales=(9, 17)):
+    gen = torch.Generator().manual_seed(seed)
+    pts = torch.rand(B, 3, generator=gen) * 2.2 - 1.1  # some points outside [-1, 1]
+    rotations = torch.rand(T, 3, 3, generator=gen) * 2 - 1
+    grids = nn.ParameterList(
+        nn.Parameter(torch.rand(3 * T, C, s, s, generator=gen) * 0.4 + 0.1) for s in scales
+    )
+    pts = pts.to(device).requires_grad_(requires_grad)
+    rotations = rotations.to(device).requires_grad_(requires_grad)
+    grids = grids.to(device)
+    if not requires_grad:
+        for g in grids:
+            g.requires_grad_(False)
+    return pts, rotations, grids
+
+
+@requires_gpu
+def test_gpu_matches_cpu_reference():
+    pts_c, rot_c, grids_c = _inputs("cpu")
+    expected = interpolate_ms_features_tilted(pts_c, grids_c, rot_c)
+    pts_g, rot_g, grids_g = _inputs("cuda")
+    actual = interpolate_ms_features_tilted(pts_g, grids_g, rot_g).cpu()
+    torch.testing.assert_close(actual, expected, rtol=1e-4, atol=1e-5)
+
+
+@requires_gpu
+def test_gpu_grads_match_cpu_reference():
+    pts_c, rot_c, grids_c = _inputs("cpu", requires_grad=True)
+    interpolate_ms_features_tilted(pts_c, grids_c, rot_c).square().sum().backward()
+    pts_g, rot_g, grids_g = _inputs("cuda", requires_grad=True)
+    interpolate_ms_features_tilted(pts_g, grids_g, rot_g).square().sum().backward()
+    torch.testing.assert_close(pts_g.grad.cpu(), pts_c.grad, rtol=1e-3, atol=1e-5)
+    torch.testing.assert_close(rot_g.grad.cpu(), rot_c.grad, rtol=1e-3, atol=1e-5)
+    for gg, gc in zip(grids_g, grids_c):
+        torch.testing.assert_close(gg.grad.cpu(), gc.grad, rtol=1e-3, atol=1e-5)
+
+
+@requires_gpu
+@requires_quantem_cuda
+def test_dispatches_to_kernel_per_scale(monkeypatch):
+    import quantem.cuda.core.ml
+
+    calls = []
+    real = quantem.cuda.core.ml.kplanes_tilted_fuse
+
+    def spy(pts, rotations, plane):
+        calls.append(plane.shape)
+        return real(pts, rotations, plane)
+
+    monkeypatch.setattr(quantem.cuda.core.ml, "kplanes_tilted_fuse", spy)
+    pts, rot, grids = _inputs("cuda")
+    interpolate_ms_features_tilted(pts, grids, rot)
+    assert len(calls) == len(grids)
+
+
+@requires_gpu
+@requires_quantem_cuda
+def test_kill_switch_forces_torch_path(monkeypatch):
+    import quantem.cuda.core.ml
+
+    def boom(pts, rotations, plane):
+        raise AssertionError("kernel should not be called with use_cuda_kernels=False")
+
+    monkeypatch.setattr(quantem.cuda.core.ml, "kplanes_tilted_fuse", boom)
+    pts, rot, grids = _inputs("cuda")
+    pts_c, rot_c, grids_c = _inputs("cpu")
+    expected = interpolate_ms_features_tilted(pts_c, grids_c, rot_c)
+    # config.set lacks __exit__, so restore explicitly rather than via `with`.
+    config.set({"use_cuda_kernels": False})
+    try:
+        actual = interpolate_ms_features_tilted(pts, grids, rot).cpu()
+    finally:
+        config.set({"use_cuda_kernels": True})
+    torch.testing.assert_close(actual, expected, rtol=1e-4, atol=1e-5)
+
+
+@requires_gpu
+def test_non_fp32_takes_torch_path():
+    pts, rot, grids = _inputs("cuda")
+    pts64 = pts.double()
+    rot64 = rot.double()
+    grids64 = nn.ParameterList(nn.Parameter(g.double()) for g in grids)
+    out = interpolate_ms_features_tilted(pts64, grids64, rot64)
+    assert out.dtype == torch.float64
diff --git a/tests/tomography/test_cuda_dispatch.py b/tests/tomography/test_cuda_dispatch.py
new file mode 100644
index 00000000..db6e9691
--- /dev/null
+++ b/tests/tomography/test_cuda_dispatch.py
@@ -0,0 +1,91 @@
+"""Tests for ``tv_loss_vol_sq`` and its optional quantem-cuda dispatch.
+
+The helper must produce the same value and gradients on every path: pure
+torch on CPU, pure torch on GPU (kill-switch or quantem-cuda absent), and
+the fused CUDA kernel when quantem-cuda is installed. The dispatch itself is
+asserted by monkeypatching the kernel entry point, so these tests are
+meaningful both with and without quantem-cuda in the environment.
+"""
+
+import pytest
+import torch
+
+from quantem.core import config
+from quantem.tomography.utils import tv_loss_vol_sq
+
+from .conftest import requires_gpu
+
+requires_quantem_cuda = pytest.mark.skipif(
+    not config.get("has_quantem_cuda"), reason="requires quantem-cuda"
+)
+
+
+def tv_vol_sq_ref(obj: torch.Tensor) -> torch.Tensor:
+    tv_d = torch.pow(obj[..., 1:, :, :] - obj[..., :-1, :, :], 2).sum()
+    tv_h = torch.pow(obj[..., :, 1:, :] - obj[..., :, :-1, :], 2).sum()
+    tv_w = torch.pow(obj[..., :, :, 1:] - obj[..., :, :, :-1], 2).sum()
+    return tv_d + tv_h + tv_w
+
+
+@pytest.mark.parametrize("shape", [(7, 6, 5), (2, 7, 6, 5)])
+def test_cpu_matches_reference(shape):
+    obj = torch.rand(shape, generator=torch.Generator().manual_seed(0))
+    torch.testing.assert_close(tv_loss_vol_sq(obj), tv_vol_sq_ref(obj))
+
+
+def test_constant_volume_is_zero():
+    assert tv_loss_vol_sq(torch.ones(8, 8, 8)).item() == 0.0
+
+
+@requires_gpu
+@pytest.mark.parametrize("shape", [(7, 6, 5), (2, 7, 6, 5)])
+def test_gpu_matches_cpu_reference(shape):
+    obj = torch.rand(shape, generator=torch.Generator().manual_seed(1))
+    expected = tv_vol_sq_ref(obj)
+    actual = tv_loss_vol_sq(obj.cuda()).cpu()
+    torch.testing.assert_close(actual, expected, rtol=1e-5, atol=1e-6)
+
+
+@requires_gpu
+def test_gpu_grad_matches_cpu_reference():
+    obj = torch.rand(6, 7, 8, generator=torch.Generator().manual_seed(2))
+    v_cpu = obj.clone().requires_grad_(True)
+    tv_vol_sq_ref(v_cpu).backward()
+    v_gpu = obj.cuda().requires_grad_(True)
+    tv_loss_vol_sq(v_gpu).backward()
+    torch.testing.assert_close(v_gpu.grad.cpu(), v_cpu.grad, rtol=1e-4, atol=1e-6)
+
+
+@requires_gpu
+@requires_quantem_cuda
+def test_dispatches_to_kernel(monkeypatch):
+    import quantem.cuda.core
+
+    calls = []
+    real = quantem.cuda.core.tv_loss_sq_3d
+
+    def spy(volume):
+        calls.append(volume.shape)
+        return real(volume)
+
+    monkeypatch.setattr(quantem.cuda.core, "tv_loss_sq_3d", spy)
+    tv_loss_vol_sq(torch.rand(4, 4, 4, device="cuda"))
+    assert len(calls) == 1
+
+
+@requires_gpu
+@requires_quantem_cuda
+def test_kill_switch_forces_torch_path(monkeypatch):
+    import quantem.cuda.core
+
+    def boom(volume):
+        raise AssertionError("kernel should not be called with use_cuda_kernels=False")
+
+    monkeypatch.setattr(quantem.cuda.core, "tv_loss_sq_3d", boom)
+    obj = torch.rand(4, 4, 4, device="cuda")
+    # config.set lacks __exit__, so restore explicitly rather than via `with`.
+    config.set({"use_cuda_kernels": False})
+    try:
+        torch.testing.assert_close(tv_loss_vol_sq(obj), tv_vol_sq_ref(obj))
+    finally:
+        config.set({"use_cuda_kernels": True})