TEST: Unit tests for splinterpolator

test/Makefile

+# 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}

test/feedsRun

+#!/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')

test/test_splinterpolator.cc

+#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-5) {
+    BOOST_CHECK_EQUAL(a.size(), b.size());
+    for (auto i = 0; i < a.size(); i++) {
+      BOOST_CHECK_SMALL(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()