diff --git a/test/Makefile b/test/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..75e2eb5ac56bd38cc1e45963f93acfcf21fb26ae
--- /dev/null
+++ b/test/Makefile
@@ -0,0 +1,42 @@
+# A Makefile for miscmaths unit tests.
+include ${FSLCONFDIR}/default.mk
+
+PROJNAME    = test-miscmaths
+TESTXFILES  = test-miscmaths
+USRINCFLAGS = -DMISCMATHS_TEST_DIRECTORY=\"$(realpath .)\"
+LIBS        = -lfsl-miscmaths -lfsl-NewNifti -lfsl-cprob
+
+# The test program can be run against
+# an in-source checkout, or against
+# an installed version of the
+# libfsl-miscmaths.so library.
+#
+# If the former, the miscmaths library
+# must have been compiled before the
+# test can be compiled.
+
+# The test program uses the Boost unit
+# testing framework, which needs librt
+# # on linux.
+SYSTYPE := $(shell uname -s)
+ifeq ($(SYSTYPE), Linux)
+LIBS  += -lrt
+RPATH := -Wl,-rpath,'$$ORIGIN/..'
+endif
+ifeq ($(SYSTYPE), Darwin)
+RPATH := -Wl,-rpath,'@executable_path/..'
+endif
+
+all: ${TESTXFILES}
+
+OBJS := $(wildcard test_*.cc)
+OBJS := $(OBJS:%.cc=%.o)
+
+# We add -I.., -L.., -Wl,-rpath so that
+# in-source builds take precedence over
+# $FSLDEVDIR/$FSLDIR
+%.o: %.cc
+	$(CXX) -I.. ${CXXFLAGS} -c -o $@ $<
+
+test-miscmaths: ${OBJS}
+	$(CXX) ${CXXFLAGS} -o $@ $^ -L.. ${RPATH} ${LDFLAGS}
diff --git a/test/feedsRun b/test/feedsRun
new file mode 100755
index 0000000000000000000000000000000000000000..3967fb951be3d2643cad094e88be59f41d070bf8
--- /dev/null
+++ b/test/feedsRun
@@ -0,0 +1,14 @@
+#!/usr/bin/env python
+
+import subprocess as sp
+import               shlex
+
+def sprun(cmd):
+    print(f'RUN {cmd}')
+    sp.run(shlex.split(cmd), check=True)
+
+
+if __name__ == '__main__':
+    sprun('make clean')
+    sprun('make')
+    sprun('./test-miscmaths --report_level=detailed')
diff --git a/test/test_splinterpolator.cc b/test/test_splinterpolator.cc
new file mode 100644
index 0000000000000000000000000000000000000000..039f138e5800c6e98015729243023b8fb7242f97
--- /dev/null
+++ b/test/test_splinterpolator.cc
@@ -0,0 +1,170 @@
+#define EXPOSE_TREACHEROUS
+
+#include "miscmaths/splinterpolator.h"
+#include "NewNifti/NewNifti.h"
+
+#include <armadillo>
+
+#include <filesystem>
+#include <stdlib.h>
+#include <vector>
+
+#define BOOST_TEST_MODULE test_splinterpolator
+#include <boost/test/included/unit_test.hpp>
+
+namespace fs  = std::filesystem;
+namespace BTF = boost::unit_test::framework;
+namespace SPL = SPLINTERPOLATOR;
+
+class TestFixture {
+
+public:
+
+  // 5x5x5 volume data to be interpolated
+  std::vector<unsigned int> dims{5, 5, 5};
+  std::vector<float>        data{
+    1,   2,   3,   4,   5,   6,   7,   8,   9,   10,
+    11,  12,  13,  14,  15,  16,  17,  18,  19,  20,
+    21,  22,  23,  24,  25,  26,  27,  28,  29,  30,
+    31,  32,  33,  34,  35,  36,  37,  38,  39,  40,
+    41,  42,  43,  44,  45,  46,  47,  48,  49,  50,
+    51,  52,  53,  54,  55,  56,  57,  58,  59,  60,
+    61,  62,  63,  64,  65,  66,  67,  68,  69,  70,
+    71,  72,  73,  74,  75,  76,  77,  78,  79,  80,
+    81,  82,  83,  84,  85,  86,  87,  88,  89,  90,
+    91,  92,  93,  94,  95,  96,  97,  98,  99,  100,
+    101, 102, 103, 104, 105, 106, 107, 108, 109, 110,
+    111, 112, 113, 114, 115, 116, 117, 118, 119, 120,
+    121, 122, 123, 124, 125
+  };
+
+  // sample points, initialised in constructor
+  size_t             npts;
+  std::vector<float> xpts;
+  std::vector<float> ypts;
+  std::vector<float> zpts;
+
+  // size of up-sampled volume, initialised in constructor
+  size_t             xsz;
+  size_t             ysz;
+  size_t             zsz;
+
+  // sampled values and gradient, initialised in constructor,
+  // populated in test cases
+  std::vector<float>              values;
+  std::vector<std::vector<float>> gradient;
+
+
+  TestFixture() {
+    // the data is up-sampled by 3x
+    for (float z = -4; z < 9.01; z+=1/3.0) {
+    for (float y = -4; y < 9.01; y+=1/3.0) {
+    for (float x = -4; x < 9.01; x+=1/3.0) {
+      if (std::abs(x) < 0.01) x = 0;
+      if (std::abs(y) < 0.01) y = 0;
+      if (std::abs(z) < 0.01) z = 0;
+      xpts.push_back(x);
+      ypts.push_back(y);
+      zpts.push_back(z);
+    }}}
+    npts = xpts.size();
+    xsz  = 40;
+    ysz  = 40;
+    zsz  = 40;
+
+    // pre-allocate space to store test outputs
+    values = std::vector<float>(npts);
+    gradient.push_back(std::vector<float>(npts));
+    gradient.push_back(std::vector<float>(npts));
+    gradient.push_back(std::vector<float>(npts));
+  };
+
+  fs::path test_data_dir() {
+    return fs::path(MISCMATHS_TEST_DIRECTORY) / "test_splinterpolator";
+  }
+
+  // load a NIfTI image from the test data directory
+  std::vector<float> load_nifti(std::string filename) {
+    char*   buf;
+    float* fbuf;
+    std::vector<NiftiIO::NiftiExtension> exts;
+
+    auto hdr = NiftiIO::loadImage(test_data_dir() / filename, buf, exts);
+
+    fbuf = reinterpret_cast<float*>(buf);
+
+    std::vector<float> values(hdr.nElements());
+
+    for (auto i = 0; i < hdr.nElements(); i++) {
+      values[i] = fbuf[i];
+    }
+
+    delete [] buf;
+    return values;
+  }
+
+  // compare two float vectors
+  void compare(std::vector<float> a, std::vector<float> b, float tol=1e-6) {
+    BOOST_CHECK_EQUAL(a.size(), b.size());
+    for (auto i = 0; i < a.size(); i++) {
+      BOOST_CHECK_CLOSE(a[i], b[i], tol);
+    }
+  }
+
+  // Run the test with a Splinterpolator instance
+  void run_test(const SPL::Splinterpolator<float>& spl) {
+
+    std::string test_name = BTF::current_test_case().p_name;
+
+    std::vector<float> deriv{0, 0, 0};
+
+    for (auto i = 0; i < npts; i++) {
+      values[i] = spl(xpts[i], ypts[i], zpts[i]);
+      spl.ValAndDerivs(xpts[i], ypts[i], zpts[i], deriv);
+      gradient[0][i] = deriv[0];
+      gradient[1][i] = deriv[1];
+      gradient[2][i] = deriv[2];
+    }
+
+    compare(values, load_nifti(test_name + "_values.nii.gz"));
+
+    // derivatives are only valid for interpolation of order > 1
+    if (spl.Order() > 1) {
+      std::vector<float> gradvals;
+      gradvals.insert(gradvals.end(), gradient[0].begin(), gradient[0].end());
+      gradvals.insert(gradvals.end(), gradient[1].begin(), gradient[1].end());
+      gradvals.insert(gradvals.end(), gradient[2].begin(), gradient[2].end());
+      compare(gradvals, load_nifti(test_name + "_gradient.nii.gz"));
+    }
+  }
+};
+
+
+BOOST_FIXTURE_TEST_SUITE(test_splinterpolator, TestFixture)
+
+
+BOOST_AUTO_TEST_CASE(order_1_3d_extrap_zeros) {
+  run_test(SPL::Splinterpolator<float>(data.data(), dims, SPL::Zeros, 1));
+}
+
+BOOST_AUTO_TEST_CASE(order_1_3d_extrap_soft_zeros) {
+  run_test(SPL::Splinterpolator<float>(data.data(), dims, SPL::SoftZeros, 1));
+}
+
+BOOST_AUTO_TEST_CASE(order_1_3d_extrap_constant) {
+  run_test(SPL::Splinterpolator<float>(data.data(), dims, SPL::Constant, 1));
+}
+
+BOOST_AUTO_TEST_CASE(order_3_3d_extrap_zeros) {
+  run_test(SPL::Splinterpolator<float>(data.data(), dims, SPL::Zeros, 3));
+}
+
+BOOST_AUTO_TEST_CASE(order_3_3d_extrap_soft_zeros) {
+  run_test(SPL::Splinterpolator<float>(data.data(), dims, SPL::SoftZeros, 3));
+}
+
+BOOST_AUTO_TEST_CASE(order_3_3d_extrap_constant) {
+  run_test(SPL::Splinterpolator<float>(data.data(), dims, SPL::Constant, 3));
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test/test_splinterpolator/data.nii.gz b/test/test_splinterpolator/data.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..1846c9d4ffaa11e52468940b18fc8eccdef90170
--- /dev/null
+++ b/test/test_splinterpolator/data.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:51ba7314b3ad94eff389b6ffa556a381809b4ff9d8b16aad87e17e516c6b026f
+size 362
diff --git a/test/test_splinterpolator/order_1_3d_extrap_constant_values.nii.gz b/test/test_splinterpolator/order_1_3d_extrap_constant_values.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..de98dc5921aee151e8a001da2ea964a6b3aa4ea4
--- /dev/null
+++ b/test/test_splinterpolator/order_1_3d_extrap_constant_values.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:335ffaaf4d9d9b1dd7bb52b8290721ee842e00983997a6b85af6df65cbdc9022
+size 4830
diff --git a/test/test_splinterpolator/order_1_3d_extrap_soft_zeros_values.nii.gz b/test/test_splinterpolator/order_1_3d_extrap_soft_zeros_values.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..80104896f05ab46c3fa0576fe05a6c19caeea62e
--- /dev/null
+++ b/test/test_splinterpolator/order_1_3d_extrap_soft_zeros_values.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4b46cde43b2b530cb3c0fa2c5cca80f9b98fde362daf87330846e0c9f6b8b1db
+size 9811
diff --git a/test/test_splinterpolator/order_1_3d_extrap_zeros_values.nii.gz b/test/test_splinterpolator/order_1_3d_extrap_zeros_values.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..15102cb5b97907d7092e37504037dd278329bf62
--- /dev/null
+++ b/test/test_splinterpolator/order_1_3d_extrap_zeros_values.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:83a4f6597266e765e8505b322dce73ffb180d9687ff35d2ed35e9a50bf67d5ca
+size 2848
diff --git a/test/test_splinterpolator/order_3_3d_extrap_constant_gradient.nii.gz b/test/test_splinterpolator/order_3_3d_extrap_constant_gradient.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..0f415cebc01d6ec3ea746c7bc2563e17645fc3b6
--- /dev/null
+++ b/test/test_splinterpolator/order_3_3d_extrap_constant_gradient.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:df6740aafbaebb0bd22298cc4e3cdd2215f1e72c4f5420bd40a1aefb38e9e6d0
+size 285334
diff --git a/test/test_splinterpolator/order_3_3d_extrap_constant_values.nii.gz b/test/test_splinterpolator/order_3_3d_extrap_constant_values.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..0f571fb783be539d274eeedf1bf56852b3b18555
--- /dev/null
+++ b/test/test_splinterpolator/order_3_3d_extrap_constant_values.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:16fe794c40b86a77a10b323d5abe1647fbe1b9ce113bf91406348f42301f0d00
+size 26849
diff --git a/test/test_splinterpolator/order_3_3d_extrap_soft_zeros_gradient.nii.gz b/test/test_splinterpolator/order_3_3d_extrap_soft_zeros_gradient.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..5cde0755ca3d8b895aea5793aa95599e70a8eaec
--- /dev/null
+++ b/test/test_splinterpolator/order_3_3d_extrap_soft_zeros_gradient.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:7df5d1f856b565d3c7fca393dee00199b79c0cf4b3bdcd73f3fdec330e23bf73
+size 138307
diff --git a/test/test_splinterpolator/order_3_3d_extrap_soft_zeros_values.nii.gz b/test/test_splinterpolator/order_3_3d_extrap_soft_zeros_values.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..4847f2dbf4acec7d6c096aa487fb40d7dc4a29cc
--- /dev/null
+++ b/test/test_splinterpolator/order_3_3d_extrap_soft_zeros_values.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:970284e1f7075083d1a1abae3a053d15ced71d4e27dcd642f756bfc7faa1c93e
+size 45702
diff --git a/test/test_splinterpolator/order_3_3d_extrap_zeros_gradient.nii.gz b/test/test_splinterpolator/order_3_3d_extrap_zeros_gradient.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..57563f1db45710b103a4d838043a688845153ef3
--- /dev/null
+++ b/test/test_splinterpolator/order_3_3d_extrap_zeros_gradient.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ff221f4112fbd85693b703fbbb01e28fc59bd638ce50ef8609b97844421398df
+size 25618
diff --git a/test/test_splinterpolator/order_3_3d_extrap_zeros_values.nii.gz b/test/test_splinterpolator/order_3_3d_extrap_zeros_values.nii.gz
new file mode 100644
index 0000000000000000000000000000000000000000..a7778757e6f777002bde5fb3f0fa887bbc49118e
--- /dev/null
+++ b/test/test_splinterpolator/order_3_3d_extrap_zeros_values.nii.gz
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4ea6adb66a31f7f7071688b14bc54b0d510447dbb3428d343f9ac7fd9f244d27
+size 8958