diff --git a/MIGRATION.md b/MIGRATION.md
index 0840873..a6cf97d 100644
--- a/MIGRATION.md
+++ b/MIGRATION.md
@@ -26,17 +26,27 @@ Each ported operation appears at three "library" levels:
 | module         | kernels | cpu-impl | cuda-impl | cpu-lib | cuda-lib | lib | CPU tested |
 |----------------|:------:|:--------:|:---------:|:-------:|:--------:|:---:|:----------:|
 | distance       |  ✓     |   ✓      |   ✓       |   ✓     |   ✓      | ✓   | **yes**    |
-| posdef         |  ✓     |   ✓      |   ✓       |   ✗     |   ✗      | ✗   | no         |
+| posdef         |  ✓     |   ✓      |   ✓       |   ✓     |   ~      | ✓   | **yes**    |
+| resize         |  ✓     |   ✓      |   ✓       |   ✓     |   ~      | ✓   | **yes**    |
+| restrict       |  ✓     |   ✓      |   ✓       |   ✓     |   ~      | ✓   | **yes**    |
+| splinc         |  ✓     |   ✓      |   ✓       |   ✓     |   ~      | ✓   | **yes**    |
 | pushpull       |  ✓     |   ✓      |   ✓       |   ✗     |   ✗      | ✗   | no         |
 | reg_field      |  ✓     |   ✓      |   ✓       |   ✗     |   ✗      | ✗   | no         |
 | reg_flow       |  ✓     |   ✓      |   ✓       |   ✗     |   ✗      | ✗   | no         |
-| resize         |  ✓     |   ✓      |   ✓       |   ✗     |   ✗      | ✗   | no         |
-| restrict       |  ✓     |   ✓      |   ✓       |   ✗     |   ✗      | ✗   | no         |
-| splinc         |  ✓     |   ✓      |   ✓       |   ✗     |   ✗      | ✗   | no         |
 | tetrahedron    |  ✓     |   ✓      |   —       |   ✗     |   ✗      | ✗   | no         |
 
-"✓" for a lib column means the dispatch layer exists; only `distance` is currently
-compiled and tested (CPU). `tetrahedron` has no cuda-impl header yet.
+"✓" for a cpu-lib/lib column means the dispatch layer exists **and is CPU-compiled+tested**.
+"~" for cuda-lib means the dispatch source is written and mirrors cpu-lib, but is **not yet
+compilable**: the corresponding `cuda-impl/<module>.h` only defines device (`CUGLOB`) kernels
+and lacks the host launchers (allocate/copy shape+stride to device, launch, forward `stream`)
+that the cuda-lib layer calls — see the "not yet fixed" section. Such modules are intentionally
+**left out of the cuda-lib Makefile `MODULES`** until the launchers exist, so the CUDA build
+keeps working. `tetrahedron` has no cuda-impl header yet.
+
+### Public API naming
+`resize`/`restrict`/`splinc` would be an illegal same-name namespace+function inside `ff::cpu`,
+so the public ops are named **`resample`** (resize), **`restriction`** (restrict), and
+**`spline_coeff`** (splinc). `restriction` accumulates into `out`, so callers pre-zero it.
 
 ## Build & test (CPU)
 
@@ -64,9 +74,37 @@ CPU path against a brute-force / reference implementation, as `test_distance.cpp
    the L1 transform was never called.
 5. **cpu-lib/Makefile, lib/Makefile** — object compile rule missing `-fPIC`, so the
    shared-library link failed (`relocation R_X86_64_PC32 … recompile with -fPIC`).
+6. **kernels/threadpool.inl** — defined non-`inline` free functions and two
+   namespace-scope globals (`internal::num_threads`, `internal::global_pool`) in a
+   header. Fine for a single-module library, but a library with ≥2 module objects
+   (e.g. distance + posdef) failed to link (`multiple definition of …`). Made the
+   helpers `inline` and wrapped the two globals in Meyers-singleton accessors (C++11
+   has no inline variables). This unblocks every multi-module library.
+7. **kernels/posdef/utils.h** — three alias templates referenced undeclared names
+   (`_as_points`, `left`/`right`, and `_return_type` with no args); any posdef include
+   failed to parse. **kernels/posdef/sym.inl** — dynamic `Sym::invert` called the
+   inherited `copy_` unqualified (two-phase lookup) → `this_type::copy_`.
+8. **cpu-impl/posdef.h** — runtime `nbatch` passed as a template arg to `index2offset`;
+   a dead static-C `else` branch that misbound the dynamic specialization under C++11
+   two-arm instantiation (guarded with `(C<0?1:C)`); `delete`→`delete[]` at 4 sites.
+9. **kernels/utils.h** — `prod<size>(x)` called `typed_prod<T,size>(x, size)` with a
+   spurious extra arg (no matching overload); triggered by `restrict::loop`.
+10. **cpu-impl/{resize,restrict,splinc}.h** — wrong include prefix `"lib/…"` →
+    `"kernels/…"`; impl namespace was plain `ff::<module>` but the kernels live in
+    `ff::cpu::` (`FF_DEVICE`) so it must be `FF_NAMESPACE_BEGIN(FF)/(FF_DEVICE)/(<module>)`
+    like distance; `index2offset_nd<ndim>()` runtime-ndim → dynamic overload;
+    `jf::has_atomic_add` → `has_atomic_add`.
 
 ## Bugs found, NOT yet fixed (need a CUDA build to verify)
 
+- **cuda-impl/{posdef,resize,restrict,splinc}.h** — provide only device (`CUGLOB`)
+  kernels and **no host launchers** analogous to distance's `CUHOST dt()`. The
+  cuda-lib layer for these modules is written but cannot compile until host wrappers
+  are added, so these modules are omitted from the cuda-lib Makefile `MODULES` for
+  now. Additionally cuda-impl uses the same wrong `"lib/…"` includes, and
+  `resize.h`'s `kernelnd` passes an undefined `x` (should be `loc`) to
+  `Multiscale<ndim>::resize`. (Tasks tracked separately.)
+
 - **cuda-impl/distance_euclidean.h** — `dt()` allocates a scratch buffer of
   `vector_size * (sizeof(offset_t) + 2*sizeof(scalar_t))` bytes, but `dt_kernel`
   addresses `z`/`d` at `buf + stride_buf * n * …` and indexes each of the
diff --git a/Makefile b/Makefile
index 1027cc3..0a2170b 100644
--- a/Makefile
+++ b/Makefile
@@ -24,7 +24,7 @@ MOSUF 	  	 = o
 SOSUF      	 = so
 SONAME     	 = soname
 OMPFLAG    	 = -fopenmp
-RPATH        = -Wl,-rpath,$$ORIGIN/../lib
+RPATH        = -Wl,-rpath,'$$ORIGIN'/../lib
 USE_OPENMP 	?= 0
 
 ########################################################################
@@ -108,7 +108,13 @@ clean-cpu:
 
 MODULES = \
 	distance \
-	posdef
+	posdef \
+	resize \
+	restrict \
+	splinc \
+	pushpull \
+	reg_field \
+	reg_flow
 
 OBJECTS  = $(addprefix $(BUILDDIR)/,$(addsuffix .$(MOSUF),$(MODULES)))
 CPPFILES = $(addsuffix .cpp,$(MODULES))
diff --git a/distance.cpp b/distance.cpp
index ebb05cc..00dad2e 100644
--- a/distance.cpp
+++ b/distance.cpp
@@ -5,7 +5,7 @@
 #include "cuda/distance.h"
 #endif
 
-using namespace FF;
+FF_NAMESPACE_BEGIN(FF)
 
 #define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
 #define IS_CPU(tensor)  (tensor.device.device_type == DLDeviceType::kDLCPU || \
@@ -123,3 +123,5 @@ void dt_mesh(
 
     throw std::invalid_argument("unsupported device");
 }
+
+FF_NAMESPACE_END(FF)
diff --git a/distance.h b/distance.h
index 6b4d6cd..f7a2dc6 100644
--- a/distance.h
+++ b/distance.h
@@ -5,6 +5,8 @@
 
 FF_NAMESPACE_BEGIN(FF)
 
+#ifndef FF_LIB_BOUND_SPLINE_T
+#define FF_LIB_BOUND_SPLINE_T
 FF_NAMESPACE_BEGIN(bound_t)
 using T = int8_t;
 static constexpr T Dynamic   = -1; ///< Used to turn-off static implementations in templated classes
@@ -30,6 +32,7 @@ static constexpr T FifthOrder    =  5;
 static constexpr T SixthOrder    =  6;
 static constexpr T SeventhOrder  =  7;
 FF_NAMESPACE_END(spline_t)
+#endif // FF_LIB_BOUND_SPLINE_T
 
 /**
  * @brief Compute the Euclidean distance transform of a tensor.
diff --git a/posdef.cpp b/posdef.cpp
index 9c64708..73d60c6 100644
--- a/posdef.cpp
+++ b/posdef.cpp
@@ -5,7 +5,7 @@
 #include "cuda/posdef.h"
 #endif
 
-using namespace FF;
+FF_NAMESPACE_BEGIN(FF)
 
 #define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
 #define IS_CPU(tensor)  (tensor.device.device_type == DLDeviceType::kDLCPU || \
@@ -136,3 +136,5 @@ void sym_invert_(
 
     throw std::invalid_argument("unsupported device");
 }
+
+FF_NAMESPACE_END(FF)
diff --git a/pushpull.cpp b/pushpull.cpp
new file mode 100644
index 0000000..4d9dd4b
--- /dev/null
+++ b/pushpull.cpp
@@ -0,0 +1,90 @@
+#include <stdexcept>
+#include "pushpull.h"
+#include "cpu/pushpull.h"
+#ifdef FF_WITH_CUDA
+#include "cuda/pushpull.h"
+#endif
+
+#define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
+#define IS_CPU(tensor)  (tensor.device.device_type == DLDeviceType::kDLCPU || \
+                         tensor.device.device_type == DLDeviceType::kDLCUDAHost)
+
+FF_NAMESPACE_BEGIN(FF)
+
+void pull(
+          DLTensor & out,
+    const DLTensor & inp,
+    const DLTensor & grid,
+          int8_t     spline,
+          int8_t     bound,
+          int8_t     extrapolate,
+          int        stream )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(out))
+        return FF_CUDA::pull(out, inp, grid, spline, bound, extrapolate, stream);
+#endif
+    if (IS_CPU(out))
+        return FF_CPU::pull(out, inp, grid, spline, bound, extrapolate, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+void push(
+          DLTensor & out,
+    const DLTensor & inp,
+    const DLTensor & grid,
+          int8_t     spline,
+          int8_t     bound,
+          int8_t     extrapolate,
+          int        stream )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(out))
+        return FF_CUDA::push(out, inp, grid, spline, bound, extrapolate, stream);
+#endif
+    if (IS_CPU(out))
+        return FF_CPU::push(out, inp, grid, spline, bound, extrapolate, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+void count(
+          DLTensor & out,
+    const DLTensor & grid,
+          int8_t     spline,
+          int8_t     bound,
+          int8_t     extrapolate,
+          int        stream )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(out))
+        return FF_CUDA::count(out, grid, spline, bound, extrapolate, stream);
+#endif
+    if (IS_CPU(out))
+        return FF_CPU::count(out, grid, spline, bound, extrapolate, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+void grad(
+          DLTensor & out,
+    const DLTensor & inp,
+    const DLTensor & grid,
+          int8_t     spline,
+          int8_t     bound,
+          int8_t     extrapolate,
+          bool       abs,
+          int        stream )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(out))
+        return FF_CUDA::grad(out, inp, grid, spline, bound, extrapolate, abs, stream);
+#endif
+    if (IS_CPU(out))
+        return FF_CPU::grad(out, inp, grid, spline, bound, extrapolate, abs, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+FF_NAMESPACE_END(FF)
diff --git a/pushpull.h b/pushpull.h
new file mode 100644
index 0000000..dc575ce
--- /dev/null
+++ b/pushpull.h
@@ -0,0 +1,118 @@
+#ifndef FF_LIB_PUSHPULL
+#define FF_LIB_PUSHPULL
+#include "dlpack.h"
+#include "defines.h"
+
+#ifndef FF_LIB_BOUND_SPLINE_T
+#define FF_LIB_BOUND_SPLINE_T
+FF_NAMESPACE_BEGIN(FF)
+
+FF_NAMESPACE_BEGIN(bound_t)
+using T = int8_t;
+static constexpr T Dynamic   = -1;
+static constexpr T Zero      =  0;
+static constexpr T Replicate =  1;
+static constexpr T DCT1      =  2;
+static constexpr T DCT2      =  3;
+static constexpr T DST1      =  4;
+static constexpr T DST2      =  5;
+static constexpr T DFT       =  6;
+static constexpr T NoCheck   =  7;
+FF_NAMESPACE_END(bound_t)
+
+FF_NAMESPACE_BEGIN(spline_t)
+using T = int8_t;
+static constexpr T Dynamic       = -1;
+static constexpr T Nearest       =  0;
+static constexpr T Linear        =  1;
+static constexpr T Quadratic     =  2;
+static constexpr T Cubic         =  3;
+static constexpr T FourthOrder   =  4;
+static constexpr T FifthOrder    =  5;
+static constexpr T SixthOrder    =  6;
+static constexpr T SeventhOrder  =  7;
+FF_NAMESPACE_END(spline_t)
+
+FF_NAMESPACE_END(FF)
+#endif // FF_LIB_BOUND_SPLINE_T
+
+FF_NAMESPACE_BEGIN(FF)
+
+/**
+ * @brief Sample ("pull") a spline-encoded volume at arbitrary coordinates.
+ *
+ * Channel-last, x-first coordinate convention:
+ *   inp  : (*batch, *inshape,  C)
+ *   grid : (*batch, *outshape, D)   with D == the spatial rank (1, 2 or 3)
+ *   out  : (*batch, *outshape, C)
+ *
+ * @param out          Output tensor (pulled samples)
+ * @param inp          Input volume (spline coefficients)
+ * @param grid         Sampling coordinates (in voxels, x-first)
+ * @param spline       Spline order applied to every spatial dim
+ * @param bound        Boundary condition applied to every spatial dim
+ * @param extrapolate  1: always; 0: not past voxel centres; -1: not past edges
+ * @param stream       Cuda stream on which to operate
+ */
+void pull(
+          DLTensor & out,
+    const DLTensor & inp,
+    const DLTensor & grid,
+          int8_t     spline      = spline_t::Quadratic,
+          int8_t     bound       = bound_t::DCT2,
+          int8_t     extrapolate = 1,
+          int        stream      = 0
+);
+
+/**
+ * @brief Splat ("push") values into a volume; numerical adjoint of `pull`.
+ *        `out` must be pre-zeroed by the caller (values are accumulated).
+ *
+ *   inp  : (*batch, *outshape, C)
+ *   grid : (*batch, *outshape, D)
+ *   out  : (*batch, *inshape,  C)
+ */
+void push(
+          DLTensor & out,
+    const DLTensor & inp,
+    const DLTensor & grid,
+          int8_t     spline      = spline_t::Quadratic,
+          int8_t     bound       = bound_t::DCT2,
+          int8_t     extrapolate = 1,
+          int        stream      = 0
+);
+
+/**
+ * @brief Splat ones (== push of an all-ones input). `out` (*batch,*inshape,1)
+ *        must be pre-zeroed.
+ */
+void count(
+          DLTensor & out,
+    const DLTensor & grid,
+          int8_t     spline      = spline_t::Quadratic,
+          int8_t     bound       = bound_t::DCT2,
+          int8_t     extrapolate = 1,
+          int        stream      = 0
+);
+
+/**
+ * @brief Sample the spatial gradients of a spline-encoded volume.
+ *
+ *   inp  : (*batch, *inshape,  C)
+ *   grid : (*batch, *outshape, D)
+ *   out  : (*batch, *outshape, C, D)
+ */
+void grad(
+          DLTensor & out,
+    const DLTensor & inp,
+    const DLTensor & grid,
+          int8_t     spline      = spline_t::Quadratic,
+          int8_t     bound       = bound_t::DCT2,
+          int8_t     extrapolate = 1,
+          bool       abs         = false,
+          int        stream      = 0
+);
+
+FF_NAMESPACE_END(FF)
+
+#endif // FF_LIB_PUSHPULL
diff --git a/reg_field.cpp b/reg_field.cpp
new file mode 100644
index 0000000..51ecc88
--- /dev/null
+++ b/reg_field.cpp
@@ -0,0 +1,55 @@
+#include <stdexcept>
+#include "reg_field.h"
+#include "cpu/reg_field.h"
+#ifdef FF_WITH_CUDA
+#include "cuda/reg_field.h"
+#endif
+
+#define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
+#define IS_CPU(tensor)  (tensor.device.device_type == DLDeviceType::kDLCPU || \
+                         tensor.device.device_type == DLDeviceType::kDLCUDAHost)
+
+FF_NAMESPACE_BEGIN(FF)
+
+void field_matvec(
+          DLTensor & out       ,
+    const DLTensor & inp       ,
+    const double   * voxel_size,
+    const double   * absolute  ,
+    const double   * membrane  ,
+    const double   * bending   ,
+          int8_t     bound     ,
+          int        ndim      ,
+          int        stream    )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(out))
+        return FF_CUDA::field_matvec(out, inp, voxel_size, absolute, membrane, bending, bound, ndim, stream);
+#endif
+    if (IS_CPU(out))
+        return FF_CPU::field_matvec(out, inp, voxel_size, absolute, membrane, bending, bound, ndim, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+void field_diag(
+          DLTensor & out       ,
+    const double   * voxel_size,
+    const double   * absolute  ,
+    const double   * membrane  ,
+    const double   * bending   ,
+          int8_t     bound     ,
+          int        ndim      ,
+          int        stream    )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(out))
+        return FF_CUDA::field_diag(out, voxel_size, absolute, membrane, bending, bound, ndim, stream);
+#endif
+    if (IS_CPU(out))
+        return FF_CPU::field_diag(out, voxel_size, absolute, membrane, bending, bound, ndim, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+FF_NAMESPACE_END(FF)
diff --git a/reg_field.h b/reg_field.h
new file mode 100644
index 0000000..16e6940
--- /dev/null
+++ b/reg_field.h
@@ -0,0 +1,90 @@
+#ifndef FF_LIB_REG_FIELD
+#define FF_LIB_REG_FIELD
+#include "dlpack.h"
+#include "defines.h"
+
+#ifndef FF_LIB_BOUND_SPLINE_T
+#define FF_LIB_BOUND_SPLINE_T
+FF_NAMESPACE_BEGIN(FF)
+
+FF_NAMESPACE_BEGIN(bound_t)
+using T = int8_t;
+static constexpr T Dynamic   = -1;
+static constexpr T Zero      =  0;
+static constexpr T Replicate =  1;
+static constexpr T DCT1      =  2;
+static constexpr T DCT2      =  3;
+static constexpr T DST1      =  4;
+static constexpr T DST2      =  5;
+static constexpr T DFT       =  6;
+static constexpr T NoCheck   =  7;
+FF_NAMESPACE_END(bound_t)
+
+FF_NAMESPACE_BEGIN(spline_t)
+using T = int8_t;
+static constexpr T Dynamic       = -1;
+static constexpr T Nearest       =  0;
+static constexpr T Linear        =  1;
+static constexpr T Quadratic     =  2;
+static constexpr T Cubic         =  3;
+static constexpr T FourthOrder   =  4;
+static constexpr T FifthOrder    =  5;
+static constexpr T SixthOrder    =  6;
+static constexpr T SeventhOrder  =  7;
+FF_NAMESPACE_END(spline_t)
+
+FF_NAMESPACE_END(FF)
+#endif // FF_LIB_BOUND_SPLINE_T
+
+FF_NAMESPACE_BEGIN(FF)
+
+/**
+ * @brief Apply a spatial regulariser operator to a multi-channel field.
+ *
+ * Layout is (*batch, *spatial, C) with `C` channels in the last axis. The
+ * penalties are per-channel weight vectors of length `C`; the highest-order
+ * non-null penalty selects the finite-difference stencil. With `voxel_size == 1`
+ * and only `absolute`, the result is a per-channel scaling
+ * `out[...,c] = absolute[c] * inp[...,c]`; with only `membrane`, it is
+ * `membrane[c]` times the discrete negative Laplacian of channel `c`.
+ *
+ * @param out         Output tensor (*batch, *spatial, C)
+ * @param inp         Input  tensor (*batch, *spatial, C)
+ * @param voxel_size  [ndim] spatial voxel size (nullptr -> all ones)
+ * @param absolute    [C] absolute (L2) penalty weights (nullptr -> zeros)
+ * @param membrane    [C] membrane penalty weights (nullptr -> disabled)
+ * @param bending     [C] bending penalty weights (nullptr -> disabled)
+ * @param bound       Boundary condition applied to every spatial dim
+ * @param ndim        Number of spatial dimensions (1, 2 or 3)
+ * @param stream      Cuda stream on which to operate
+ */
+void field_matvec(
+          DLTensor & out       ,
+    const DLTensor & inp       ,
+    const double   * voxel_size = nullptr,
+    const double   * absolute  = nullptr,
+    const double   * membrane  = nullptr,
+    const double   * bending   = nullptr,
+          int8_t     bound     = bound_t::DCT2,
+          int        ndim      = 1,
+          int        stream    = 0
+);
+
+/**
+ * @brief Diagonal (preconditioner) of the regulariser operator, same
+ *        conventions as `field_matvec`. Writes into `out` (*batch, *spatial, C).
+ */
+void field_diag(
+          DLTensor & out       ,
+    const double   * voxel_size = nullptr,
+    const double   * absolute  = nullptr,
+    const double   * membrane  = nullptr,
+    const double   * bending   = nullptr,
+          int8_t     bound     = bound_t::DCT2,
+          int        ndim      = 1,
+          int        stream    = 0
+);
+
+FF_NAMESPACE_END(FF)
+
+#endif // FF_LIB_REG_FIELD
diff --git a/reg_flow.cpp b/reg_flow.cpp
new file mode 100644
index 0000000..9846fea
--- /dev/null
+++ b/reg_flow.cpp
@@ -0,0 +1,55 @@
+#include <stdexcept>
+#include "reg_flow.h"
+#include "cpu/reg_flow.h"
+#ifdef FF_WITH_CUDA
+#include "cuda/reg_flow.h"
+#endif
+
+#define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
+#define IS_CPU(tensor)  (tensor.device.device_type == DLDeviceType::kDLCPU || \
+                         tensor.device.device_type == DLDeviceType::kDLCUDAHost)
+
+FF_NAMESPACE_BEGIN(FF)
+
+void flow_matvec(
+          DLTensor & out       ,
+    const DLTensor & inp       ,
+    const double   * voxel_size,
+          double     absolute  ,
+          double     membrane  ,
+          double     bending   ,
+          int8_t     bound     ,
+          int        ndim      ,
+          int        stream    )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(out))
+        return FF_CUDA::flow_matvec(out, inp, voxel_size, absolute, membrane, bending, bound, ndim, stream);
+#endif
+    if (IS_CPU(out))
+        return FF_CPU::flow_matvec(out, inp, voxel_size, absolute, membrane, bending, bound, ndim, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+void flow_diag(
+          DLTensor & out       ,
+    const double   * voxel_size,
+          double     absolute  ,
+          double     membrane  ,
+          double     bending   ,
+          int8_t     bound     ,
+          int        ndim      ,
+          int        stream    )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(out))
+        return FF_CUDA::flow_diag(out, voxel_size, absolute, membrane, bending, bound, ndim, stream);
+#endif
+    if (IS_CPU(out))
+        return FF_CPU::flow_diag(out, voxel_size, absolute, membrane, bending, bound, ndim, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+FF_NAMESPACE_END(FF)
diff --git a/reg_flow.h b/reg_flow.h
new file mode 100644
index 0000000..0d37852
--- /dev/null
+++ b/reg_flow.h
@@ -0,0 +1,90 @@
+#ifndef FF_LIB_REG_FLOW
+#define FF_LIB_REG_FLOW
+#include "dlpack.h"
+#include "defines.h"
+
+#ifndef FF_LIB_BOUND_SPLINE_T
+#define FF_LIB_BOUND_SPLINE_T
+FF_NAMESPACE_BEGIN(FF)
+
+FF_NAMESPACE_BEGIN(bound_t)
+using T = int8_t;
+static constexpr T Dynamic   = -1;
+static constexpr T Zero      =  0;
+static constexpr T Replicate =  1;
+static constexpr T DCT1      =  2;
+static constexpr T DCT2      =  3;
+static constexpr T DST1      =  4;
+static constexpr T DST2      =  5;
+static constexpr T DFT       =  6;
+static constexpr T NoCheck   =  7;
+FF_NAMESPACE_END(bound_t)
+
+FF_NAMESPACE_BEGIN(spline_t)
+using T = int8_t;
+static constexpr T Dynamic       = -1;
+static constexpr T Nearest       =  0;
+static constexpr T Linear        =  1;
+static constexpr T Quadratic     =  2;
+static constexpr T Cubic         =  3;
+static constexpr T FourthOrder   =  4;
+static constexpr T FifthOrder    =  5;
+static constexpr T SixthOrder    =  6;
+static constexpr T SeventhOrder  =  7;
+FF_NAMESPACE_END(spline_t)
+
+FF_NAMESPACE_END(FF)
+#endif // FF_LIB_BOUND_SPLINE_T
+
+FF_NAMESPACE_BEGIN(FF)
+
+/**
+ * @brief Apply a spatial regulariser operator to a vector flow field.
+ *
+ * Layout is (*batch, *spatial, C) with `C == ndim` flow components in the last
+ * axis. The operator is the sum of the requested penalties (absolute, membrane,
+ * bending); the highest-order non-zero penalty selects the finite-difference
+ * stencil. With `voxel_size == 1` and only `absolute`, the result is
+ * `absolute * inp`; with only `membrane`, it is `membrane` times the discrete
+ * negative Laplacian of the field.
+ *
+ * @param out         Output tensor (*batch, *spatial, ndim)
+ * @param inp         Input  tensor (*batch, *spatial, ndim)
+ * @param voxel_size  [ndim] spatial voxel size (nullptr -> all ones)
+ * @param absolute    Absolute (L2) penalty weight
+ * @param membrane    Membrane (first-order) penalty weight
+ * @param bending     Bending (second-order) penalty weight
+ * @param bound       Boundary condition applied to every spatial dim
+ * @param ndim        Number of spatial dimensions (1, 2 or 3)
+ * @param stream      Cuda stream on which to operate
+ */
+void flow_matvec(
+          DLTensor & out       ,
+    const DLTensor & inp       ,
+    const double   * voxel_size = nullptr,
+          double     absolute  = 0.0,
+          double     membrane  = 0.0,
+          double     bending   = 0.0,
+          int8_t     bound     = bound_t::DCT2,
+          int        ndim      = 1,
+          int        stream    = 0
+);
+
+/**
+ * @brief Diagonal (preconditioner) of the regulariser operator, same
+ *        conventions as `flow_matvec`. Writes into `out` (*batch, *spatial, ndim).
+ */
+void flow_diag(
+          DLTensor & out       ,
+    const double   * voxel_size = nullptr,
+          double     absolute  = 0.0,
+          double     membrane  = 0.0,
+          double     bending   = 0.0,
+          int8_t     bound     = bound_t::DCT2,
+          int        ndim      = 1,
+          int        stream    = 0
+);
+
+FF_NAMESPACE_END(FF)
+
+#endif // FF_LIB_REG_FLOW
diff --git a/resize.cpp b/resize.cpp
new file mode 100644
index 0000000..18aeff0
--- /dev/null
+++ b/resize.cpp
@@ -0,0 +1,34 @@
+#include <stdexcept>
+#include "resize.h"
+#include "cpu/resize.h"
+#ifdef FF_WITH_CUDA
+#include "cuda/resize.h"
+#endif
+
+#define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
+#define IS_CPU(tensor)  (tensor.device.device_type == DLDeviceType::kDLCPU || \
+                         tensor.device.device_type == DLDeviceType::kDLCUDAHost)
+
+FF_NAMESPACE_BEGIN(FF)
+
+void resample(
+          DLTensor & out    ,
+    const DLTensor & inp    ,
+          int8_t     spline ,
+          int8_t     bound  ,
+          double     shift  ,
+    const double   * scale  ,
+          int        ndim   ,
+          int        stream )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(out))
+        return FF_CUDA::resample(out, inp, spline, bound, shift, scale, ndim, stream);
+#endif
+    if (IS_CPU(out))
+        return FF_CPU::resample(out, inp, spline, bound, shift, scale, ndim, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+FF_NAMESPACE_END(FF)
diff --git a/resize.h b/resize.h
new file mode 100644
index 0000000..09cb365
--- /dev/null
+++ b/resize.h
@@ -0,0 +1,71 @@
+#ifndef FF_LIB_RESIZE
+#define FF_LIB_RESIZE
+#include "dlpack.h"
+#include "defines.h"
+
+#ifndef FF_LIB_BOUND_SPLINE_T
+#define FF_LIB_BOUND_SPLINE_T
+FF_NAMESPACE_BEGIN(FF)
+
+FF_NAMESPACE_BEGIN(bound_t)
+using T = int8_t;
+static constexpr T Dynamic   = -1;
+static constexpr T Zero      =  0;
+static constexpr T Replicate =  1;
+static constexpr T DCT1      =  2;
+static constexpr T DCT2      =  3;
+static constexpr T DST1      =  4;
+static constexpr T DST2      =  5;
+static constexpr T DFT       =  6;
+static constexpr T NoCheck   =  7;
+FF_NAMESPACE_END(bound_t)
+
+FF_NAMESPACE_BEGIN(spline_t)
+using T = int8_t;
+static constexpr T Dynamic       = -1;
+static constexpr T Nearest       =  0;
+static constexpr T Linear        =  1;
+static constexpr T Quadratic     =  2;
+static constexpr T Cubic         =  3;
+static constexpr T FourthOrder   =  4;
+static constexpr T FifthOrder    =  5;
+static constexpr T SixthOrder    =  6;
+static constexpr T SeventhOrder  =  7;
+FF_NAMESPACE_END(spline_t)
+
+FF_NAMESPACE_END(FF)
+#endif // FF_LIB_BOUND_SPLINE_T
+
+FF_NAMESPACE_BEGIN(FF)
+
+/**
+ * @brief Resample (prolongation) a tensor to a new shape using spline
+ *        interpolation. Prefilter the input with `spline_coeff` first for a
+ *        proper interpolating resize.
+ *
+ * The `ndim` trailing dimensions are spatial; leading dimensions are batch.
+ * For each output voxel `loc`, input is sampled at `x = scale*(loc+shift)-shift`.
+ *
+ * @param out     Output tensor (*batch, *outshape)
+ * @param inp     Input tensor  (*batch, *inshape)
+ * @param spline  Spline order applied to every spatial dim
+ * @param bound   Boundary condition applied to every spatial dim
+ * @param shift   Anchor shift (0 centres, 0.5 edges)
+ * @param scale   [ndim] per-dim scaling (input-index per output-index)
+ * @param ndim    Number of spatial dimensions (1, 2 or 3)
+ * @param stream  Cuda stream on which to operate
+ */
+void resample(
+          DLTensor & out    ,
+    const DLTensor & inp    ,
+          int8_t     spline = spline_t::Quadratic,
+          int8_t     bound  = bound_t::DCT2,
+          double     shift  = 0.0,
+    const double   * scale  = nullptr,
+          int        ndim   = 1,
+          int        stream = 0
+);
+
+FF_NAMESPACE_END(FF)
+
+#endif // FF_LIB_RESIZE
diff --git a/restrict.cpp b/restrict.cpp
new file mode 100644
index 0000000..6681ce5
--- /dev/null
+++ b/restrict.cpp
@@ -0,0 +1,34 @@
+#include <stdexcept>
+#include "restrict.h"
+#include "cpu/restrict.h"
+#ifdef FF_WITH_CUDA
+#include "cuda/restrict.h"
+#endif
+
+#define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
+#define IS_CPU(tensor)  (tensor.device.device_type == DLDeviceType::kDLCPU || \
+                         tensor.device.device_type == DLDeviceType::kDLCUDAHost)
+
+FF_NAMESPACE_BEGIN(FF)
+
+void restriction(
+          DLTensor & out    ,
+    const DLTensor & inp    ,
+          int8_t     spline ,
+          int8_t     bound  ,
+          double     shift  ,
+    const double   * scale  ,
+          int        ndim   ,
+          int        stream )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(out))
+        return FF_CUDA::restriction(out, inp, spline, bound, shift, scale, ndim, stream);
+#endif
+    if (IS_CPU(out))
+        return FF_CPU::restriction(out, inp, spline, bound, shift, scale, ndim, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+FF_NAMESPACE_END(FF)
diff --git a/restrict.h b/restrict.h
new file mode 100644
index 0000000..26d3674
--- /dev/null
+++ b/restrict.h
@@ -0,0 +1,70 @@
+#ifndef FF_LIB_RESTRICT
+#define FF_LIB_RESTRICT
+#include "dlpack.h"
+#include "defines.h"
+
+#ifndef FF_LIB_BOUND_SPLINE_T
+#define FF_LIB_BOUND_SPLINE_T
+FF_NAMESPACE_BEGIN(FF)
+
+FF_NAMESPACE_BEGIN(bound_t)
+using T = int8_t;
+static constexpr T Dynamic   = -1;
+static constexpr T Zero      =  0;
+static constexpr T Replicate =  1;
+static constexpr T DCT1      =  2;
+static constexpr T DCT2      =  3;
+static constexpr T DST1      =  4;
+static constexpr T DST2      =  5;
+static constexpr T DFT       =  6;
+static constexpr T NoCheck   =  7;
+FF_NAMESPACE_END(bound_t)
+
+FF_NAMESPACE_BEGIN(spline_t)
+using T = int8_t;
+static constexpr T Dynamic       = -1;
+static constexpr T Nearest       =  0;
+static constexpr T Linear        =  1;
+static constexpr T Quadratic     =  2;
+static constexpr T Cubic         =  3;
+static constexpr T FourthOrder   =  4;
+static constexpr T FifthOrder    =  5;
+static constexpr T SixthOrder    =  6;
+static constexpr T SeventhOrder  =  7;
+FF_NAMESPACE_END(spline_t)
+
+FF_NAMESPACE_END(FF)
+#endif // FF_LIB_BOUND_SPLINE_T
+
+FF_NAMESPACE_BEGIN(FF)
+
+/**
+ * @brief Restriction: the adjoint (transpose) of the resize prolongation.
+ *
+ * Maps a fine tensor `inp` to a coarse tensor `out`; `out` is ACCUMULATED into
+ * and must be zeroed by the caller. Use the reciprocal scale of the matching
+ * resize (resize: coarse/fine; restriction: fine/coarse) with the same shift.
+ *
+ * @param out     Output (coarse) tensor (*batch, *outshape), pre-zeroed
+ * @param inp     Input (fine) tensor    (*batch, *inshape)
+ * @param spline  Spline order applied to every spatial dim
+ * @param bound   Boundary condition applied to every spatial dim
+ * @param shift   Anchor shift (0 centres, 0.5 edges)
+ * @param scale   [ndim] per-dim scaling (input-index per output-index)
+ * @param ndim    Number of spatial dimensions (1, 2 or 3)
+ * @param stream  Cuda stream on which to operate
+ */
+void restriction(
+          DLTensor & out    ,
+    const DLTensor & inp    ,
+          int8_t     spline = spline_t::Quadratic,
+          int8_t     bound  = bound_t::DCT2,
+          double     shift  = 0.0,
+    const double   * scale  = nullptr,
+          int        ndim   = 1,
+          int        stream = 0
+);
+
+FF_NAMESPACE_END(FF)
+
+#endif // FF_LIB_RESTRICT
diff --git a/splinc.cpp b/splinc.cpp
new file mode 100644
index 0000000..f0df084
--- /dev/null
+++ b/splinc.cpp
@@ -0,0 +1,30 @@
+#include <stdexcept>
+#include "splinc.h"
+#include "cpu/splinc.h"
+#ifdef FF_WITH_CUDA
+#include "cuda/splinc.h"
+#endif
+
+#define IS_CUDA(tensor) (tensor.device.device_type == DLDeviceType::kDLCUDA)
+#define IS_CPU(tensor)  (tensor.device.device_type == DLDeviceType::kDLCPU || \
+                         tensor.device.device_type == DLDeviceType::kDLCUDAHost)
+
+FF_NAMESPACE_BEGIN(FF)
+
+void spline_coeff(
+          DLTensor & inp_out ,
+          int8_t     spline  ,
+          int8_t     bound   ,
+          int        stream  )
+{
+#ifdef FF_WITH_CUDA
+    if (IS_CUDA(inp_out))
+        return FF_CUDA::spline_coeff(inp_out, spline, bound, stream);
+#endif
+    if (IS_CPU(inp_out))
+        return FF_CPU::spline_coeff(inp_out, spline, bound, stream);
+
+    throw std::invalid_argument("unsupported device");
+}
+
+FF_NAMESPACE_END(FF)
diff --git a/splinc.h b/splinc.h
new file mode 100644
index 0000000..29fb042
--- /dev/null
+++ b/splinc.h
@@ -0,0 +1,62 @@
+#ifndef FF_LIB_SPLINC
+#define FF_LIB_SPLINC
+#include "dlpack.h"
+#include "defines.h"
+
+#ifndef FF_LIB_BOUND_SPLINE_T
+#define FF_LIB_BOUND_SPLINE_T
+FF_NAMESPACE_BEGIN(FF)
+
+FF_NAMESPACE_BEGIN(bound_t)
+using T = int8_t;
+static constexpr T Dynamic   = -1;
+static constexpr T Zero      =  0;
+static constexpr T Replicate =  1;
+static constexpr T DCT1      =  2;
+static constexpr T DCT2      =  3;
+static constexpr T DST1      =  4;
+static constexpr T DST2      =  5;
+static constexpr T DFT       =  6;
+static constexpr T NoCheck   =  7;
+FF_NAMESPACE_END(bound_t)
+
+FF_NAMESPACE_BEGIN(spline_t)
+using T = int8_t;
+static constexpr T Dynamic       = -1;
+static constexpr T Nearest       =  0;
+static constexpr T Linear        =  1;
+static constexpr T Quadratic     =  2;
+static constexpr T Cubic         =  3;
+static constexpr T FourthOrder   =  4;
+static constexpr T FifthOrder    =  5;
+static constexpr T SixthOrder    =  6;
+static constexpr T SeventhOrder  =  7;
+FF_NAMESPACE_END(spline_t)
+
+FF_NAMESPACE_END(FF)
+#endif // FF_LIB_BOUND_SPLINE_T
+
+FF_NAMESPACE_BEGIN(FF)
+
+/**
+ * @brief In-place spline coefficient prefilter along the last dimension.
+ *
+ * Filters the last axis (all leading axes are batch); prefiltering makes
+ * spline interpolation reproduce the input samples. To filter several axes,
+ * permute the tensor and call repeatedly.
+ *
+ * @param inp_out  Input/Output tensor in DLTensor format (float32/float64)
+ * @param spline   Spline order (orders 0/1 are no-ops)
+ * @param bound    Boundary condition
+ * @param stream   Cuda stream on which to operate
+ */
+void spline_coeff(
+          DLTensor & inp_out ,
+          int8_t     spline   = spline_t::Cubic,
+          int8_t     bound    = bound_t::DCT2,
+          int        stream   = 0
+);
+
+FF_NAMESPACE_END(FF)
+
+#endif // FF_LIB_SPLINC